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Preface 


Tackling this thesis project required a search for knowledge among several 
areas previously beyond my grasp. These included neurophysiology, mathematics, 
pattern recognition and artificial neural network culture, and software development. 
Thanks to Dr. Steven K. Rogers, Dr. Matthew Kabrisky, Captain Dennis Ruck and 
Captain Greg Tarr, I had access to expertise in all of these disciplines. I would never 
have conceived of this research project, let alone attempted a crack at it, without 
their guidance. I am further indebted to my wife, Tonda, for her support end 
understanding throughout my pursuit of the degree requirements. This experience 
has immeasurably improved my knowledge of science and engineering, also leaving 
me with an appreciation for many complex concepts beyond my comprehension. 
As expected, earning this MSEE cdso entailed plenty of grief and aggravation. In 
fact, I almost turned down this rare educational opportunity before it began, until I 
pondered its monetary cost. It was free. 


Laurence Edward Lazofson 
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Abstract 

This thesis project included a literature survey of biological and artificial neu¬ 
ral network research followed by development and testing of high-order and image 
recognition hierarchical neural network algorithms. Following training, performance 
testing of second-order and third-order networks yielded maximum accuracies com¬ 
parable to those achieved by multilayer perceptron classifiers operating on test data 
sets. Several versions of an image classification algorithm were tested for learning 
performance using pixel data from forward-looking infrared (FLIR) images of tanks, 
trucks, target boards, <ind clutter. Employing the biologically-motivated Lamberti- 
zation and contrast normalization of pixel windows, correlations with multiple Gabor 
function wavelets, and a “phase synchronizing” local averaging routine, the image 
classification network extracted data features. Different network versions fed the 
extracted features to varying output classification schemes. To improve separation 
of problem classes, recommendations were made for varying the parameters of the 
Gabor function wavelets and modifying the phase synchronization scheme to extract 
more suitable features from image pixel data. 
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A BIOLOGICALLY-INSPIRED 
NEURAL NETWORK ARCHITECTURE 
FOR IMAGE PROCESSING 


L Introduction 

1.1 Summary of Current Knowledge 

It is widely known that in the striate (visual) cortex, certain key features of 
visual image data are segmented following mapping from the retinas. These sepa¬ 
rately distinguished aspects of visual information include intensity (of three colors), 
motion, binocular disparity, and texture (14). Texture is essentially characterized 
by spatial frequency, a term often used to describe both pitch and orientation of 
edges. In other words, texture consists of locdized sinusoidal grating patterns, or in¬ 
tensity functions, with varying periodicity and angular orientation in the 2-D visual 
field (13). Figure 1 illustrates this concept. 

Seeking a better understanding of cortical processes, neurophysiological re¬ 
searchers have investigated the functioning of corticeil tissue. As a “crumpled sheet” 
of neural mass, the mammalian cortex is divided into cortical columns. This 2-D 
array of adjacent cortical columns functionally resembles a honeycomb (13). In the 
visual cortex (Area 17), each cortical column is selectively tuned to respond to a 
different aspect of segmented visual data. Blasdel and Sadama employed “voltage- 
sensitive dyes” to observe “detailed maps of orientation selectivity,” or functional 
organization of cortical columns, in the visuail cortex of macaque monkeys (2). 

Experimental observations of neuronal functioning in the visual cortex aue used 
in this thesis to develop a biologically-motivated artificial neural network system for 
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Figure 1. 2-D Windowed Patterns of Varying Periodicities and Orientations 

image processing. Observed neuronal processes highlighted in this reseeirch included 
stimulus-specific responsiveness (4, 9, 10, 11, 12, 26), phase synchronous firing of 
groups of neurons (4, 9, 26), the semblance of Gabor functions in neuronal receptive 
field plots (4, 11, 12), and “axo-axonic interconnections” (19) in biological neural 
networks. The artificial neural network system developed in this thesis emulated 
the characteristics of these observed biological phenomena. Details eure discussed in 
Chapter II. 

1.1.1 Neuronal Stimulus Specificity Ervin, Gray and Singer, Hubei and 
Wiesel, and Jones, Stepnoski, and Palmer independently observed stimulus-specific 
responsiveness of neurons in cat visual cortex (4, 9, 10, 12). Each research team 
varied stimuli exposed to the visual field of test subjects and measured neuron firing 
outputs. One unified result from these researchers seems to verify that the visual 
cortex is an org 2 inized collection of columns, with some of these columns sensitive to 
specific orientations of a stimulus (14). 
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1.1.2 Phase Synchrony Gray and Singer further observed synchronous oscil¬ 
lations of stimulus-specific neurons in the visual cortex of cats (9). These findings 
indicate a possible cortical mechanism that enables neurons responding to unique 
stimuli to simultaneously broadcast their message to higher cortical levels (26). Such 
a mechanism facilitates the “binding” of local stimulus details into larger, global fea¬ 
ture aspects. Attempting to mathematically model these findings, Kammen, Holmes, 
and Koch proposed and tested two algorithms for phase locking artificial cortical 
colimms (15). These concepts are discussed in greater detail in Chapter II of this 
thesis. 

1.1. S Gabor Functions Modeling Neuronal Responses Several biological re¬ 
searchers also recorded neuronal receptive field responses that mapped to the gen¬ 
eral form of Gabor functions (Figure 2) (4, 11, 12). This may be a highly significcint 
finding, as a Gabor elementary function might model a transfer function of a neuron 
being stimulated within a lineaur region of operation (14). 

An elementary Gabor function consists of a sinusoidal function multiplied (or 
windowed) by a Gaussian function. A Gabor function is consequently chajracter- 
ized by the frequency of the sinusoidal cofactor and the variance of the Gaussian 
cofactor. The Gaussian function localizes, or truncates, the sinusoid. The sinusoid 
discriminates frequency. (1) 

Independent research by Ayer and Fretheim implemented 2-D spatial Gabor 
functions in segmenting image features. Both researchers correlated image scenes 
with multiple Gabor functions of differing spatial frequencies and variemces (1, 6). It 
is possible that the visual cortex employs a similar correlation process in discerning 
texture information on an image (14). 

1 . 1.4 “Axo-Axonic” Neuronal Connections In addition to the well-known 
“axo-dendritic” neural connection model depicted in Figure 3, neurophysiological 
resecLTch also indicates the existence of “axo-aixonic” synapses in biological neural 
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Figure 2. Examples of 2-D Gabor Functions (1:3) 


networks (19:2). These connections suggest processing interactions, or modulations, 
among neuron outputs. The “axo-axonic” connection serves as a model for high- 
order artificial neural networks discussed later in this document (19). 

1.1.5 Theory The pursuit of this thesis work required tmdem consideration 
of previously developed artificial neural network algorithms and the cited biolog¬ 
ical observations of neural functioning in cat visual cortex. In the visual cortex, 
phase synchronous firing of stimulus-specific groups of neurons may indeed tempo¬ 
rally unite their conveyed messages. This would enable simultaneous summation 
and processing of these dynamic neuronal responses at higher cortical levels (26). 
However, in implementing a software model of a network, computer memory storage 
essentially freezes time by holding calculated numerical outputs of artificial neurons. 
A subsequent step then sums these stored values in other Mtificial neurons (typ¬ 
ically in a higher layer). Software modeling thus removes the temporal aspect of 
neuronal oscillatory firing by storing numerical output values that represent firing 









late. This; skirts the need for synchronizing outputs of software neurons, whereas 
phase locking phenomena seem crucial in biological neural networks with transient 
oscillatory outputs. Similarly, it is questionable whether biologically measured tem¬ 
poral latency, an aspect of phase delay, would prove valuable in software models of 
neural networks (23). This concept is discussed further in Chapter II. 

Results from cited reseaichers also allude to neuronal stimulus specificity that 
function 2 Jly responds in the semblance of imiqwe Glabor spatial mappings. It is 
widely believed that neurons in the visual cortex operate as Gabor detectors, re¬ 
sponding to unique combinations of pitch and orientation to segment localized image 
textures (14). This implication influenced the direction of this thesis work. 

1.2 Background 

1.2.1 Artificial Neural Networks For several decades, researchers have been 
developing and testing a variety of artificial neural network algorithms for use on 
certain classes of problems. These algorithms, whether supervised, unsupervised, or 
hybrid models, are capable of “learning” functional relationships among data sets. 
Network learning, also referred to as “training,” occurs as the network is presented 
with examples of vector feature data. For some cleisses of networks, the training 
process drives iterative mathematical adjustments of node connection “weights,” 
or functional multipliers, until network convergence. Following training, a viable 
network should be able to “generalize,” or solve a problem, with new, unknown data 
inputs. (23) 

Unlike conventional electronic digital computer serial processing, neural net¬ 
work algorithms employ densely interconnected nodes modeling massive parallelism. 
As seen in biological neural networks, massively parallel architectures yield a more 
robust system, characterized by graceful degrzwlation rather than being prone to to¬ 
tal failure due to oue m.alfunctioning computational element (18, 23). Ultimately, 
software-developed neural networks may be holographically based or implemented 
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electronically on silicon chips. 

Emulating some of the observed euspects of biological neurons (Figure 3) and 
neural networks, artificial neural networks have been designed to solve classifica¬ 
tion, identification ^md optimization problems, perform functional estimation and 
prediction, echeince signal-to-noise ratios, and serve in a variety of other applica¬ 
tions (3, 23). Image and speech recognition and prediction of chaotic time series are 
specific examples of potential applications of artificial neural networks. 

1.2.2 High-Order Artificial Neural Networks High-order neural network al- 
goritiiins have been shown to surpass the performance of “classicad” multilayer feed¬ 
forward networks for some applications (8,19, 20). By feeding polynomial functional 
combinations of vector data components into a single layer oi artificial neurons, these 
networks may efficiently “carve out” complex, functional decision regions in multi¬ 
dimensional feature spaces. High-order networks can thus tackle problems that are 
not linearly separable. Compared to multilayer first-order networks, tests of these 
high-order models indicated faster convergence and greater probability of attaining 
solutions for some problems (8, 19, 20). 

Multilayer first-order networks (Figure 4) also functionally adapt to nonlinear¬ 
ities by implicitly implementing higher orders in their hidden node layers. However, 
the added network layers «md additional weight modifications may yield slower, and 
less probable, convergence to solutions. For these reasons, high-order network al¬ 
gorithms have been pursued in solving more complex, nonlinearly-separable prob¬ 
lems. (8, 21) 

l.S Problem Statement 

There existed a need to perform a consolidated literature review of biological 
and artificial neural network research. Incorporating the concepts and phenomena 
culled from this integrated literature survey, there w?s a subsequent need to expand 
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Figure 3. Emulating the Biological Neuron (5:8) 









Figure 4. Example of Multilayer Network of Artificial Neurons (23) 









the state of the art of existing neural network algorithms. 


1.4 Objective 

This thesis work sought to survey, integrate, and document a subset of bio¬ 
logical and artificial neural network research. Having initially reconciled this body 
of knowledge, this thesis subsequently aimed to develop and test a software model 
consistent with the cited research findings. This entailed emulating some of the bi¬ 
ological observations of cortical functioning. The intent was to employ high-order 
networks and implement Gabor functions zmd the concept of phase synchrony to 
create a better model for solving certain classes of engineering problems, including 
image processing. The model ultimately employed Lambertization and contrast nor¬ 
malization of windowed image regions, correlations with multiple Gabor functions, a 
phase synchronizing local averaging routine, and a feedforward network cleissification 
layer. Details are covered in Chapter III. 

1.5 Scope and Limitations 

This work did not attempt to provide a comprehensive, unified explanation 
of cortical processing. Rather, this thesis intended to describe amd emulate some 
of the observed characteristics of neuronal functioning in the visuad cortex, amd 
extend previously developed artificial neural network models. This work did not 
address the validity of biological observations and proposed cortical phenomena by 
cited researchers. This thesis did not assess whether some of the described cortical 
measurements may have been artifacts. 

This reseau-ch focused on a subset of observed neural characteristics, specifically 
citing work by Ervin (4), Gray and Singer (9), and Jones, Stepnoski, and Palmer (11, 
12). Attempting to unify the biological observations of these researchers, this thesis 
subsequently sought to develop an artificial neur 2 d network model consistent with 
the cited aspects of neurophysiology. This model was related to current network 
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algorithms and applied to the target recognition problem for forwaxd-looking infrared 

(FLIR) images. 

1.6 Definitions 

1- D One-dimensional. 

2- D Two-dimensional. 

Artifact An observed measurement not related to a process or phenomena of inter¬ 
est, typically resulting from a random or systematic error in experimental or 
measurement procedures. 

Artificial Neural Network (Net) Software or haxdwaure emulation of simple, ob¬ 
servable characteristics of biological neural networks. 

Artificial Neuron A computational element that mathematically models a biolog¬ 
ical neuron. 

Chaotic Time Series A deterministic function that is seemingly random (23). 

Gabor Function A sinusoidal fimction multiplied (or windowed by) a Gaussian 
function, characterized by the frequency of the sinusoidal cofactor and the 
variance of the Gaussian cofactor (1). 

Neuron Nerve cell. 

Node A single computational element of a neural network. 

Receptive (Visual) Field The visual region over which a specific neuron responds 
(in the visual cortex) (10:109). 

Segment To distinguish features or aspects of interest from background. 

Spatial FVequency Pitch (periodicity) and angular orientation of a function or 
structure of interest, such as an image. 

Striate Cortex Section of cortex that receives mapped visual information, also 
known as visual cortex, VI, or Area 17 (14). 
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Weights Multiplicative factors betweeu nodes or between inputs and nodes, auedo- 
gous to the strength of synaptic connections between neurons (17). Also called 
“connection weights.” 

Wetware Biological cell computational elements (23). 

1.7 Sequence of Prtsentation 

This introductory chapter discussed the problem and background concepts for 
this thesis. Chapter II provides the literature review for this research. Chapter III 
presents the research approach, discussing methods and algorithms. Results and 
recommendations are covered in Chapter IV amd Chapter V, respectively. 
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IL Literature Review 


Several of the published research findings cited in this thesis depict various 
experimental observations of neuronal processes in the striate (visual) cortex. These 
observations of appjirent cortical functioning include stimulus-specific responsive¬ 
ness and phase synchrony of neuronal firing, the semblance of Gabor functions in 
characterizing neuronal receptive fields, and “axo-axocic interconnections” (19) in 
biologicsd neural networks. A key aspect of this thesis work, addressed in this chap¬ 
ter, aimed to reconcile the different cited experimental observations. The subsequent 
research goal sought the development of biologically-motivated neural network algo¬ 
rithms for solving certain classes of engineering problems. For this thesis work, the 
problem of interest was image processing. 

8.1 Neurophysiological Research 

8 .1.1 Ervin Research published by Ervin in 1965 (4) provides a solid basis 
for stimulus specificity of neuronal response, phase synchronicity in groups of firing 
neurons, and use of Gabor functions in chzuracterizing neuronal receptive field pro¬ 
files. Ervin’s paper helps unify these aspects of cortical processing into a coherent 
view supported by more current research articles. For this reason, his relatively 
early discoveries are discussed along with recently published research. Ervin’s ex¬ 
perimented observations followed Hubei and Wiesel’s widely acknowledged research 
characterizing neuronal receptive fields and firing responses in cat visual cortex (10). 

Using a computer display to generate a variety of visual stimuli m front of 
the eyeball of an immobilized cat, Ervin experimentally tracked response outputs in 
the cat’s visual cortex. This was done with a microelectrode and a gross electrode, 
both inserted into the visual cortex and positioned relatively close together. Limited 
in size, the microelectrode measured firing output of a single neuron responding to 
each visual stimulus. T1 j larger gross electrode, piercing through many neurons in 
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the surrounding cortical region, measured the average evoked response to the same 
visual stimulus. Data were stored by computer. (4) 

Unfortunately, Ervin’s article did not explicitly state whether or not the cat 
was anesthetized. In assessing the validity of measured neuronal responses, this 
may be a significant factor, as neurons in an anesthetized subject might not display 
typical behavior (13). 

While monitoring the firing rate of a single neuron with the microelectrode, 
Ervin initially measured the neuron’s output response as a function of position in 
its visual receptive field, A small spot of light, analogous to an impulse input, 
was shifted point by point across the visual field. The neuronal firing rate was 
measured for each field position of the coristant-intensity spot of light (4). Ervin 
employed this initial procedure to locate the “visual field center” (4:37), pinpointing 
the location in the neuron’s receptive field where the light spot induced the greatest 
firing response. Over twenty years later, Jones and Palmer performed a similar 
experimental procedure to characterize the spatial receptive field response of neurons 
in cat striate cortex (11). 

Having located the visual field center, Ervin employed it as a center point for 
visual stimuli. Applied stimuli included lines of varying angular orientation, circles, 
squares, and triangles of different sizes, and multiple point patterns. Each unique 
stimulus Bashed before the cat’s eye for a duration of 400 ms. The microelectrode and 
gross electrode recorded the firing responses for each test neuron and its rurrounding 
region, respectively. (4:36-38) 

For each specific stimulus, Ervin plotted the individual neuron’s firing rare, 
the average evoked response of the airea surrounding the neuron, and the stimulus 
duration, all versus time (Figure 5). Actually, he divided the time axis into discrete 
20 ms interval bins. For each interval, the number of firings of the test neuron were 
counted and plotted (4:39-40). The resulting histogram effectively conveys firing 
frequency versus time, though in an indirect way. 
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Figure 5. Neuronal Firing Histogram (Top), Stimulus Duration (Middle), and Av¬ 
erage Evoked Response (Bottom) Plots with 20ms Bins (4:40) 


Comparisons of the resulting one-dimensional plots revealed a startling discov¬ 
ery. A myriad of different plot patterns for varying stimuli indicated stimulus-specific 
responses in neurons (4, 9,10,11,12, 26). Each stimulus generated a unique pattern 
of “phasic,” “tonic,” and inhibitory components along the time scale (4:40). As used 
by Ervin, phasic components referred to relatively transient peaking points on the 
firing rate plots, while tonic components described plateau regions (4). Ervin’s ob¬ 
servations of stimulus specific responsiveness followed previous work by Hubei and 
Wiesel, who claimed “shape, position, and orientation” (10:151) of stimuli affect 
neuronaJ firing responses. 

Ervin conveyed another crucial concept by juxtaposing the plotted response of 
each neuron with the average evoked response of its surrounding region. In genercd, 
the plot patterns for both the individued neuron and the average evoked response 
were similar (4:40) Figure 5 conveys this idea. This may indicate synchronous 
oscillatory activity of groups of neurons (4, 9, 26). Such an inference, if true, may 
be of tremendous significance in understanding higher cortical processes. Recent 
research by Gray and Singer further suggested this phenomena (9). 


14 





Figure 6. l-D Plots of Neuronal Receptive Field Responses (4:43) 

Similar to his procedure for initially determining a neuron’s visual field center, 
Ervin collected additional data to spatially chairacterize a neuron’s entire receptive 
field for different phasic components of its firing response. This process essentially 
froze time at the occurrence of a phasic component on a neuron’s plot of firing rate 
versus time. Firing rate was then plotted as a function of spatial position across 
the neuron’s 2-D visual receptive field (4:42-44). Ervin illustrated this relationship 
on both 1-D and 2-D spatial coordinate plots (Figure 6). It is interesting to note 
that these plots resemble the envelopes of 1-D Gabor functions. This finding is in 
accord with later work by Jones amd Palmer (11). Ervin observed that different 
phasic components (etirly and late) resulting from a single stimulus produced unique 
receptive field response mappings (4:42). 










Ervin’s plots of firing rate versus time jdso indicated that responses to simul¬ 
taneous stimuli at multiple points in a neuron’s receptive field did not necessarily 
sum (4:45). This suggests that, lacking the property of superposition, a neuron is 
not obliged to respond linearly to stimuli. This should preclude unrestricted use of a 
previously characterized impulse response to predict cell, or system, output (13, 21). 

Depending on a given neuron’s “choice” of specific stimulus, however, it might 
linearly add the correct combination of inputs (13, 21). For example, a neuron with 
a “preference” for edges at a given orientation may respond linearly to two points 
of light extrapolating to a line following this “preferred” orientation. Otherwise, 
if two light points form a segment of different orientation in the visual field, this 
same neuron may be inhibited and behave nonlinearly. Ervin suggests this idea of 
orientational preference (4:50). Furthermore, the linearity assumption for neural 
response may be valid for limited ranges of stimuli (14). 

S.1.8 Jones and Palmer In 1987, Jones and Palmer employed a procedure 
they termed “reverse correlation” (11) to characterize the 2-D spatial receptive field 
responses of neurons in cat striate cortex. Similar to Ervin's results, their clear, 
2-D mappings of the visual receptive field response closely approximated the form 
of Gabor functions. 

Fourteen anesthetized, adult cats were exposed to visual stimuh on an oscillo¬ 
scope screen. The stimuli consisted of small, rectangular, bright and dark spots pre¬ 
sented for varying time durations, typically 50 or 100 ms (11). The reverse correla¬ 
tion process electroniceilly stored, and later matched, mea.sured neuron spike outputs 
with the appropriate stored stimuli generated on the oscilloscope screen (11:1191). 
Receptive field functions for bright and dark stimuli were then plotted separately 
(firing response versus 8pati^Ll field position). To fuse the receptive field plots for 
both bright and dairk stimuli, Jones and Palmer subtracted the dark stimulus out¬ 
put functions from corresponding light stimulus output functions. Thus the dark 
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Figure 7. 2-D Spatial Plots of Neuronal Receptive Field Responses (11:1196) 

stimulus output patterns emulated inhibition of brightness in the 2-D spatial do¬ 
main (11:1192). As noted by Jones and Palmer (11), the resulting receptive field 
response functions appeared to be Gabor functions (Figure 7). 

As Jones and Pcdmer indicated, previous work by other researchers attempting 
to map neuronal receptive fields with one-dimensional measurements might prove in¬ 
valid if the overall response functions were not Cartesian separable (11:1188). Rather 
than independently sweep across length and width axes, Jones and Palmer charac¬ 
terized the receptive fields in two spatial dimensions. This method provided accurate 
2-D representations of the field mappings, regeirdless of iheir functional Cartesian 
separability (11:1188). In fact, Jones and Palmer ascertained that approximately 
half of the neurons investigated displayed 2-D receptive field functions that were not 
Cartesian separable (11:1208). 
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Figure 8. Varying Delay Times for 2-D Spatial Plots of Neuronal Receptive Field 
Response (11:1197) 

Jones and Palmer further employed the reverse correlation process to assess 
the time delay betwreen stimuli and corresponding output spikes in the test neuron. 
They varied the delay time at values of 0, 60, 100, and 150 ms in the reverse corre¬ 
lation algorithm. Their graphical results showed that at a time delay of 50 ms after 
stimulus presentation, a distinguishable functional form of the neuronal receptive 
field response emerged (Figure 8). At time delays of 0, 100, Jind 150 ms after stim¬ 
ulus presentation, the receptive field plots depicted mostly noise with no apparent 
response function. (11:1195-1197) 

Jones’ and Palmer’s observations can be reconciled with Ervin’s findings. For 
a given test neuron and a specific applied stimulus, emergence of a receptive field 
response, in the form of a Gabor function, occurred 50 ms after stimulus presenta¬ 
tion (11:1195-1197). This pomt in time marked what Ervin called a phasic compo¬ 
nent, or transient r,iaxima of neuronal firing, indicating the time until appearance 
of the receptive field function. In light of Ervin’s findings, the unique receptive 
field response profile and its 50 ms delay are both probably unique to the stimulus 
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pieseEted and the test neuron (4, 11). 

It would be interesting to see an extension of Jones’ and Palmer’s results 
through continued use of their reverse correlation process beyond a 150 ms delay 
time. Jones and Palmer observed the emergence of a receptive field response, in 
the form of a Gabor function, 50 ms following the visual stim?ilus. They discontin¬ 
ued the reverse correlation process after a 150 ms delay (11:1197). For the applied 
stimulus, there may yet have existed additional, later phasic components for the 
neuronal response, as Ervin observed on firing rate versus time plots (4:40). These 
later components, perhaps emerging at delays of 200 ms and beyond, may give rise 
to different receptive field response mappings, although still perhaps in the general 
form of Gabor functions. It is possible that different combinations of test neurons 
and stimuli might yield unique delay times for one or more phasic components, each 
with unique 2-D Gabor field mappings. Extending the reverse correlation proce¬ 
dure in this manner could confirm Ervin’s observations and expand understanding 
of stimulus specificity in neuronal responses. 

S.l.S Jones, Stepnoski, and Palmer Beyond the previously described re¬ 
search characterizing 2-D spatial receptive field responses of neurons, Jones, Step- 
noski, and Palmer examined the 2-D spectral response of 36 neurons in cat visual 
cortex (12:1212). Their experimental procedure was similar to that of the Jones and 
Palmer experiment cited earlier. In this experiment, however, the visual stimuli were 
sinusoidal gratings drifting across the neuronal receptive field (12:1212). 

Varying the spatial frequency and orientation of these 2-D sinusoidal inten¬ 
sity functions, Jones et al meeisured and graphicedly depicted neuronal response 
(spiking) as a fimction of spatial frequency and orientation appearing in the visual 
receptive field (12). Jones et al noted that a corticed cell’s response to the drifting 
gratings yielded “a rectified sinusoidal modulation of the spike frequency” where 
‘'the degree of rectification v£iried from cell to cell, but for each cell, the form of 
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Figure 9. 2-D Spectral Plot of Neuronal Receptive Field Response (12:1222) 

the response was constant irrespective of stimulus spatiad frequency, orientation, or 
contrast” (12:1212). 

Examining the plots in this article, these neuronal spectral responses appeared 
in the form of Fourier transforms of Gabor functions, manifested as two, separate 2- D 
Gaussian functions (Figure 9). This is consistent, since Gabor functional forms in the 
spatial receptive field should correspond to Fourier transforms of these functions in 
the spectral domain. Mathematically described, a sinusoid multiplied by a Gaussian 
function Fourier transforms to two impulses in the frequency domain convolved with 
a Gaussian (22). This convolution yields two Gaussians. Figure 10 illustrates this 
mathematical process. It should also be noted that the spectred response plots were 
polar and were obtained assuming polar separability of the orientation ((?) and spatial 
frequency (p) one-dimensional functions (12:1223). 

These results further support the idea of stimulus specificity of neurons in 
the visual cortex (4, 9, 10, 11, 12, 26). Responses to varying sinusoidal grating 
stimuli indicate that neurons segment on the image features of spatial frequency and 
orientation (12). Jones et al appeirently define spatial frequency as pitch distance 
between edges. In this definition, they do not include the angular orientation of the 
parallel grating edges, considering orientation a separate feature (12). 
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Figure 10. 2-D Plots Depicting Fourier TVansforming of a Gabor Function; F Mod 

els Neuronal Spectral Response (6) 
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8.1.4 Gray and Singer In 1988, Gray and Singer documented experimen¬ 
tal data suggesting that “local neuronal populations in the visual cortex engage in 
stimulus-specific synchronous oscillations” (9:1698). Their findings, in accord with 
other research cited here, were obtained from cortical measurements on 15 adult 
cats and 12 kittens, all anesthetized (9:1698). Stimuli exposed to the visual recep¬ 
tive fields consisted of illuminated line segments at varying orientations, velocities, 
and direction of movement (9:1699). 

Gray and Singer observed stimulus-specific responses in the test neurons (9). 
Perhaps even more impressive was their observation that, for a given stimulus, “ad¬ 
jacent neurons” fired “simultaneously and in synchrony” (9:1701). Equally signif¬ 
icant was their observation that different groups of firing neurons, “spatially sep¬ 
arated” (9:1702) on the cortex by up to 7 mm, oscillated in phase, provided they 
selected for the “same orientation specificity” (26:298). 

Assessing Gray and Singer’s results, Stryker summarized, “receptive fields of 
the neurons at the two sites had a common orientation specificity and were aligned 
so that they could be stimulated by a single long bar of light” (26:298). Stryker 
further suggested that the “global property of the stimulus” (26:298) affected firing 
correlations between spatially separated groups of neurons. He concluded this since 
one long light bar, rather tham two separate light bars of the saime orientation, yielded 
synchronous firings of spatially separated cortical areas (26;298). Evidently, the two 
separate bars “did not bridge the gap between the two receptive fields” (26:298). 

Gray and Singer suggested that “adjacent neurons” firing “simultaneously and 
in synchrony” for a specific stimulus are “confined approximately to a single orienta¬ 
tion column” (9:1701). They also stated that, based on measurements with multiple 
electrodes, “synchronization across spatially separate columns does occur” (9:1702). 
Stryker deduced that many “neurons in the visual cortex,” responding simultane¬ 
ously to the same stimulus, can better convey their message to higher cortical areas 
by broadcasting “in unison” (26:297). Such phenomena may provide ideas for an 
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enhanced neural network model as “the phase of the oscillatory response may be 
used as a further dimension of coding” (9:1702). The concept of phase synchrony 
essentially “binds” locad stimulus details into lau'ger, global feature aspects. 

2.2 Biologically-Motivated Neural Network Models 

2.2.1 Kammen, Holmes, Knch Based on Gray amd Singer’s observations of 
corticad processes, Kaunmen, Holmes, amd Koch subsequently proposed two algo¬ 
rithmic models to mathematically simulate phaise locking among groups of neu¬ 
rons (15). Their first model incorporated a “one-dimensional array” of symbolic 
cortical columns, each “coupled” to its two “neairest neighbors” (15:1-182). Their 
second model connected each cortical column to “a common compaurator which feeds 
back a function of the average phase” (15:1-182). See Figure 11 for these represen¬ 
tations. 

Anadyzing their first model, the one emulating “nearest neighbor couplings,” 
Kaimmen et al determined that it did not successfully phase lock mathematical 
inputs (cortical columns) “separated by laurge numbers of inactive (unstimulated) 
units” (15:1-182). Increaising the number of corticad units, or lengthening the one- 
dimensional auray, cam be expected to exacerbate the problem. It is apparent that 
this algorithm did not accurately model the biological phenomena observed by Gray 
amd Singer (15). 

The second model, however, successfully coupled the phaise and firing frequency 
of the excited (stimulated) cortical columns, amd moved them out of phase (and off 
frequency) with the unstimulated columns (15:1-182). More succinctly, Kammen et al 
determined that “not only do the excited units fire at the same rate, but they remaun 
exactly in phase regardless of their geometrical arrangement or separation” (15:T 
182). Recall that this model implemented a common comparator for average phaise 
feedback (15:1-182). These results offer further direction for a phase-synchronous, 
enhamced haudware neural network model. 
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Figure 11. Phase Synchronous Models: Ne^u'est Neighbor and Common Compara¬ 
tor (15:1-183) 
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2.2.2 Daugman Recognizing the importance of Gabor functions in modeling 
neuronal reci^ptive field responses, Daugm 2 m developed a three-layer neural network 
to output Gabor transform coefficients for input images. The image inputs were 2-D 
intensity functions mapped across an array of pixels. (3) 

According to Daugman, typical images are characterized by localized regular¬ 
ities in texture, brightness, and edge continuation. These consistencies make such 
images good candidates for minimizing their coding information. Daugman’s neural 
network accomplishes this by transforming a 2-D pixel array intensity function into 
Gabor coefficients. This algorithm reduces the amount of electronic data storage 
needed to describe and reproduce an image. (3) 

2.2.S High-Order Artificial Neural Networks 

2.2.3.1 Basic Theory High-order networks feed vector inputs via high- 
order polynomial functions into a single layer of artificial neurons. Compared to 
standard multilayer networks, high-order networks more efficiently “carve out” com¬ 
plex, functional decision regions within feature spaces for some problems. Stand£ird 
first-order backpropagation networks implicitly implement higher orders in their liid- 
den layers of artificial neurons. However, the standard feedforward networks’ added 
layers and additional weight modifications may result in slower convergence toward 
problem solutions. (8, 21) 

2.2.5.2 Meador Assuming the existence of axon-to-axon connections 
in biological neural networks, Meador developed single-layer, high-order, artificial 
neural networks to model modulations, oi multiplicative interactions, among ax¬ 
ons. Meador indicated that axo-axonic interconnections are physiologically docu¬ 
mented, yet have not been widely discussed in artificial neural network literature. 
For relative simplicity, Meador limited the scope of his models to implement the 
widely recognized “axo-dendritic” and the lesser acknowledged “axo-axo-dendritic 
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synapses.” (19) 


Limiting his mathematical algorithms, Meador modeled second-order connec¬ 
tions in his networks, although higher order connections are possible to imple¬ 
ment (28). Meador’s first algorithm, an axo-axonic model with zuljustable weights 
for linear and second-order polynomial input functions, appeeired as the following 
equation (19:3): 


I I 

neti = 53 di,(l + 53 aijkXk)xj (1) 

In Meador’s equation shown above, 1 < * < where I is the number of inputs, 
N is the number of neurons, each r is an input vector component value, d,j the linear 
connection weights, ^uld the higher order (axo-axonic) weight factors. Meador 
revised this equation to an equivalent form he termed “a second-order quadratic 
interconnect” (19:4): 


I II 

H DijXj + 53 X] ^iikXkXj (2) 

J=1 J=1 k=l 

In comparison to the first equational model, this equation is characterized by 
a different weight space where Dij = dij and Aijk = di^a^k (19:4). To modify the 
network weights, Meador used a gradient-descent method. This mathematical rule 
minimized the squared error (difference) between the desired output and the actual 
output for input vectors. 

Testing the high-order networks on three quadiatically separable problems (in¬ 
cluding the XOR problem), Meador found that the axo-axonic weight set converged 
to solutions faster than the model implementing quadratic interconnect weights. 
However, both of these single-layer, high-order networks converged faster, and with 
greater probability, than a standard, multilayer backpropagation network with no 
high-order inputs (19:15). 
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2.2.S.S Tenorio and Lee Similar to Meador’s model, Tenorio and Lee 
developed a “Self Organizing Neural Network” (28). This algorithm modified in¬ 
ternal weights, and selected appropriate linear £ind higher order functions to solve 
identification and classification problems (28:57). They tested the algorithm with a 
MacKay-Glaiss differential equation (28). 

Tenorio and Lee did not use a mean-square error method for weight adjust¬ 
ment. Instead, supervised network learning was achieved with a “modified Minimum 
Description Length (MDL) criterion,” also referred to as “Structure Estimation Cri¬ 
terion (SEC)” (28:58). This weight modification scheme mathematically selected a 
network function ba«''d an simplicity and best estimate. The transfer function of 
the network wa.« c^Tnposed of a combination of line2U‘ and higher order function2il 
inputs. (28:58) 

As the mathematical basis for their algorithm, Tenorio and Lee cited the 
Kolmogorov-Garbor polynomial, as shown in the following equation (28:57): 

y = ao + X] + Z] Z) + - (3) 

• ♦ i 

Variables x and y represent system input and output, respectively (28:58). Note the 
mathematical similarity to Meador’s previously cited approach (19). Like Meador, 
Tenorio and Lee opted to simplify their model by limiting high-order functions to 
quadratic polynomials (28:58). 

To update weights, the model was trained with the following equation for 
MDL (28:60): 

MDL = - log /(x I ^) 4- 0.5A: log N (4) 

“where /(x | 0) is the estimated probability density function of the model, k is the 
number of parameters, and N is the number of observations” (28:60). Tenorio and 
Lee contended this was a viable model based on successful prediction of the MacKay- 
Gliiss equation (a chaotic time series). They ^dso noted the trade-off between level 
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of complexity «uid modeling accuracy for the edgorithm (28:63). 


2 . 2 .5 .4 Zhang and Miller Zhang and Miller proposed a high-order net¬ 
work scheme aimed at improving a model developed in 1987 by Heilenberg. Heilen- 
berg’s model was designed “to explain how sensory maps could enhance resolution 
through orderly arraingement of broadly tuned receptors” (29:444). In essence, the 
model was intended to mathematically emulate the phenomenon of hyperacuity, a 
condition attained by a system resolving inputs to detciil finer than “inter-receptor 
spacing” (29: 444 - 44 5). 

While Heilenberg’s model employed only lineeu: weights, Zhang and Miller in¬ 
cluded Hermitian polynomial weighting functions, converting the Heilenberg algo¬ 
rithm to a higher order version (29:445-446). They showed that the improved algo¬ 
rithm successfully modeled hyperacuity for an eirray of receptor inputs with a defined 
inter-receptor spacing (29). 

The general form of the Hermitian polynomial, where p is the order of 

the polynomial, is represented as: 

= (5) 

2.2.3.5 Namatame Naunatame presented a unique high-order network 
employing Chebychev polynomials as inputs to a multilayer system (Figure 12). The 
network was designed to learn nonlinear continuous functions. Namatame contended 
that “first-order multi-layer networks and the high-order flat networks without hid¬ 
den units ... are inadequate to generadize amd to learn the nonlinear structures un¬ 
derlying the continuous mappings” (20:1-682). Namatame demonstrated the use of 
nonmonotonic Chebychev polynomials feeding three layers of aurtificiaJ neurons (20:1- 
682). 

Namatame is in accord with previously cited research regarding the improved 
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Figure 12. A Multilayer High-Order Neural Network (20:1-683) 


capabilities of high-order networks over standard first-order, multilayer models. How¬ 
ever, Namatame’s use of multiple layers in conjunction with high-order functional 
inputs deviates from the “classical” single slab versions of higb-order implemen¬ 
tations. As the results suggest, for learning and predicting nonlinear continuous 
functions, the added aspect of multiple layers may prove a superior approach for 
high-order networks. (20) 

The Chebychev polynomials (Figure 13), of order i, aie defined for 0 < x < 1 
in the following general equation (20:1-681): 

, , cos (iarccos(2i — 1))-f-1 . , 

9ii^) = -^^- - - ( 6 ) 
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figure 13. Examples of Chebychev Polynomials 
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\ 2.2.3.6 Le Cun An image recognition algorithm, tested ou handwritten 

lumerical digits, was developed by Le Cun et al. The algorithm employed multiple 
layers of nodes in a network scheme designed for “extracting local features and 
combining them to form higher-order features” (16:44). Handwritten digits for the 
data set were obtained from zip code numereds off addressed envelopes processed by 
the US Postal Service. The training set was comprised of 7291 digit exemplars while 
2007 digits were segregated for the test set (16:41). 

Each input image was mapped on a 16-by-16 pixel array. The first layer of 
nodes was organized into twelve groups of 8-by-8 node arrays. Weights were con¬ 
strained for each group of nodes, with each node “viewing” a 5-by-5 pixel block from 
the input image. With 64 constrained weights per node group, and unconstrained 
thresholds, weight updating was performed via backpropagation using a computed 
average error (16). This architecture is similar to Fukushima’s hiercirchical neocog- 
uitron except that it is supervised and not binary (7). 

The next layer of nodes consisted of twelve groups of 4-by-4 node arrays, each 
node being fed from a 5-by-5 node block from the previous layer. Outputs from this 
layer were fully connected to a layer of thirty hidden nodes, each subsequently feeding 
all ten output nodes (16). Figure 14 diagrams the network architecture developed 
by Le Cun et al. Following treiining, this network achieved 95% accuracy on the test 
data (16:45). This algorithm also inspired, in part, tls proposed image recognition 
network developed in this thesis and detailed in Chapter III. 

2.3 Literature Review Summary 

Neurophysiological reseeurchers cited in this literature reviev/ have observed a 
variety of function 2 d measurements in the neomammalian visual cortex. Their obser¬ 
vations lend credence to the proposed phenomena of stimulus specificity and phase 
synchronous firing of groups of neurons, and the potential use of Gabor functions 
in modeling neuronal (spatial) receptive field profiles. Meador, whose high-order 
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Figure 14. Le Cun’s Digit Recognition Network (16:44) 
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networks outperformed standard first-order networks, also cited biological literature 
supporting axon-to-axon connections among neurons. It is possible, however, that 
some of these biological measurements are artifacts representing some random or 
systematic noise in the experimental processes, rather than observations of actual 
cortical phenomena. 

To segr rent image features, Ayer and Fretheim independently correlated image 
scenes with Gabor functions of differing spatial frequencies and vaxiances (1, 6). The 
visu 2 d cortex may possibly employ a similar process for discerning texture (14). 

Based on the assumption that the cited biological observations represented 
real processes in the “wetware,” this thesis research sought to explore new aspects of 
neural network algorithms with ideas gained from measured cortical processes. This 
did not aissume a comprehensive understanding of cortical functioning. Rather, this 
work attempted to unify a limited body of biological, mathematical, and network 
algorithmic knowledge to develop an improved model for solving certain classes of 
engineering problems. 

This chapter discussed and integrated experimental work published by several 
biological and axtifici 2 d neural, network researchers. The following chapter covers the 
methodology and approach used in this thesis. 
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III. Methodology 


Citing the work of several neurophysiologiced researchers, Chapter II covered 
a subset of observed biological processes in the mammalian visual cortex along with 
recent research into artificial neural network algorithms. Specifically noted were 
the phenomena of stimulus-specific responsiveness and phase-synchronous firing of 
neurons in the visujid cortex, the semblance of Gabor functions in characterizing 
neuronal receptive fields, and “axo-axonic interconnections” (19) in biological neurcJ 
networks. This thesis work sought to emulate some of the simple aspects of these 
biological phenomena in softwatfe algorithms designed for feature extraction and 
image classification (target identification). 

This chapter covers the development and testing of proposed software-based 
algorithms. One of the discussed solution methodologies employs multiple Gabor 
function correlations across image pixel arrays. This process, intended for feature 
extraction within an image, may emulate stimulus-specific responsiveness of some 
cortical neurons to unique local spatial frequencies. Also discussed is a local av¬ 
eraging routine intended to model phase-synchronous neuronal firing in the visual 
cortex. 

The first section of this chapter outlines the design and testing of high-order 
classification engines. Operating on extracted features of segmented images, these 
classifiers are theoretically analogous in function to higher cortical associative pro¬ 
cesses employing axo-axonic connections among neurons. 

S.l Investigation of High-order Neural Network Classifiers 

The first phase of this work entailed the development and testing of single¬ 
layer, high-order artificial neural networks (see Appendix A for software code). This 
involved the software implementation of second and third-order network algorithms. 
Written in standard C programming language and compiled and run on two ELXSI 
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Figure 15. Two Segregated Classes of the 2-D XOR Problem 

6400 computers, these networks were tested on three separate classification problems 
to assess their learning performance (speed of solution convergence and maximum 
achieved accuracies). The purpose of this segment of the research was to determine 
the potential usefulness of high-order neural networks as classifiers for a variety of 
problems, including image recognition. 

The test problems consisted of data sets for the 2-D exclusive-or (XOR) prob¬ 
lem (Figure 15), the 2-D “mesh” problem (Figure 16), and Ruck’s calculated mo¬ 
ments of pixel data from military vehicle images. The mesh data, consisting of 1000 
vectors, was obtained from Tarr, whose previous work at AFIT used the data to test 
multilayer perceptroa classifiers. The image moment data, consisting of 81 vectors, 
was computed by Ruck from pixel data and used in testing the capabilities of mul¬ 
tilayer perceptrou classifiers by Ruck Jmd Tarr at AFIT. The computed moments 
were position, scale, and rotation invariant. 
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3.1.1 Algorithm for the High-order Classifiers Tbe networks developed for 
the second and third-order classifiers were similar with a few exceptions. The second- 
order network fed linear and second-order combinations of feature inputs, each miilti- 
plied by separate weighting factors, to the sigmoid function comprising each artificial 
neuron. The third-order network fed linear, second-order, and third-order combina¬ 
tions of feature inputs to each artificial neuron. The Kolmogorov-G^lrbor polynomial, 
shown in Equation 3, depicts a generalized mathematicid version of the summation of 
these high-order input combinations. The polynomial summation, added to a thresh¬ 
old vaJue (0) for each artificial neuron, was fed to the sigmoid function internal to 
each artificial neuron. 

Second and third-order combinations of feature inputs Me simply the multiplied 
exponential and cross-product combinations of the vector components. The sigmoid 
functions to which the summation of these first, second, cind third-order combinations 
were fed may be represented by: 

1 

J'" " 1 + c-(o+«n) 

where g„ is the sigmoid output and 0n the threshold for the nth artificial neuron. 
The variable a represents a truncated form of the Kolmogorov-Garbor polynomial. 
For the second-order network algorithm: 

• » J 

For the third-order network algorithm: 

a = E + E £ H WijkXiXjXk (9) 

I i j % j k 

The array variable x represents the vector components, or feature inputs, and u? is a 
corresponding array of weight factors for feature inputs and second and third-order 


( 8 ) 


(7) 
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multiplied combinations of the feature inputs. 

The software programs were written allowing the user to specify key parame¬ 
ters for each compiled run. These variable p 2 urameters included the number of output 
classes for the problem, the number of exemplar vectors, training vectors, and com¬ 
ponents (features) per vector as dictated by the data set, the number of training 
iterations run between each accuracy test of the network, and the total number of 
training iterations to be run and plotted. Since these networks consisted of a single 
layer of artificisd neurons, the number of output classes, number of output neurons, 
and total number of neurons were all equal. 

The user also specified the learning rate, 77 , cind the name of the data file 
to be accessed by the program. There was also a “switch” in the program which 
could be set to divide the learning rate by the fan in, or number of inputs fed to 
each artificial neuron’s sigmoid function. Dividing rf by the fan in may aid network 
learning when relatively large numbers of inputs are fed to the artificicd neurons ( 21 ). 
When not selected, the fan in vjuriable waus bard set to one so that tj was not divided 
by the number of inputs. This was the case for test running the 2-D XOR 2 uid mesh 
problems. For the twenty-two-dimensional Ruck data, however, rj was divided by the 
fan in due to the vastly increased number of sigmoid inputs created by the higher 
dimensionality of the data vectors. 

After opening the data file and loading into storage airrays the flag (index) 
number, feature component V 2 jues, and desired output class for each vector, a pro¬ 
gram module calculated the number of second-order 2 uad third-order multiplicative 
combinations for each vector’s components. These cedculated numbers, used for dy¬ 
namically allocating storage 2 u:rays for second and third-order input combinations, 
are fimctions of the number of features per vector. Having correctly allocated the 
array sizes, the program multiplied out the second and third-order combinations of 
the features and stored these high-order input values in the arrays. 
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Figure 17. Second-order Network Algorithm with Two Feature Inputs 

To save memory, the algorithm avoids the commutative redundancy of multi¬ 
plied input combinations. For example, for two feature components, xj and X 2 , the 
second-order network stores and uses the inputs Xi, X 2 , xf, x^, and Xi*X 2 (Figure 17). 
Each squared exponential is calculated and used only once and the multiplied combi¬ 
nation X 2 *xi is not calculated. By eliminating the redundemcy of commutativity, the 
algorithm saves a huge number of floating-point operations emd substantial memory, 
particularly when operating on data vectors of high dimensionality. 

Initial testing of the network algorithms on the mesh problem data resulted 
in erratic fluctuations of accuracy during the leau’ning process. This occurred us¬ 
ing learning rates of 0.35 and 1.0. The seemingly random vacillations iu learned 
accuracy indicated that perhaps the networks were not learning to any degree. It 
was postulated that the disappointing results were due to data inputs that were 
not norm 2 ilized, including the higher-order multiplied combinations derived from the 
flrst-order inputs. Second and third-order multiplications of fractioneJ input values 
yield relatively small high-order inputs. Miniscule high-order inputs may contribute 
little to weight updating and learning, thus the need for normalization. 
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Normalization program modules were subsequently developed and inserted into 
the networks to verticedly normalize the first, second, and third-order vector data 
arrays. The mean, ft, and standard deviation, <t, were calculated for each feature 
component column of the exemplar vectors. All first, second, and third-order vector 
components, x,-, for exemplar and test vectors, were then vertically normalized using 
the equation (21): 

OC ■ ““ IX 

Xi(normaliz€d) = — - (10) 

(X 

After completing the initialization routine of the high-order algorithms, the 
array of 6 values, one for each neuron, and the arrays of weight values, one for each 
first, second, 2 uid third-order input for each neuron, were filled with random floating 
point numbers between -0.5 and 0.5. These vaJues lie within the region of greatest 
change in the sigmoid function. For each network run, the random function generator 
was itself randomly seeded using the computer’s time clock. This ensured a new set 
of initial weights and thresholds for the start of each learning run. 

During the training loops of the second and third-order algorithms, exemplar 
vectors from the data set were selected in random sequence for exposure to the 
network. For each exemplar vector shown to the network, inputs consisting of the 
vector’s components and their higher-order combinations were multiplied by corre¬ 
sponding weight factors and fed to the sigmoid functions comprising each neuron. 
Based on the input vector, this yielded a unique output value, j/, for each neuron. A 
desired output, d, was assigned to each neuron such that d = 1 only for the neuron 
corresponding to the output class of the presently run exemplar vector. For the 
other neurons in the layer, each of which selected for other cleisses, the algorithm set 
d = 0. For these single-layer networks, the nimiber of neurons equaled the number 
of output classes for a problem. A neuron would output high to indicate network 
selection of its corresponding output class. 

The error, Atn, was cadculated for each neuron using the actual neuron output, 
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y, and its desired output, d, for each exemplar vector. The well-known gradient error 
equation for the sigmoid function is; 

Aw = y*(l-y)*{d—y) (11) 

For each neuron, each weight value, w,, was updated via the equation: 

Wi{updated) = t«,- -}-»/ ♦ Aw * i, (12) 

where x, represents the corresponding input for the weight (including higher-order 
inputs for updating higher-order weights). As previously mentioned, rj may be di¬ 
vided by the fan in to aid network learning, particul£U‘ly when many inputs au:e fed 
to the artificial neurons (21), Threshold values (9) for each neuron were updated 
using the same equation, where the value of x,- was set to one. 

Through successive training iterations, or sequenti£il random exposures to ex¬ 
emplar vectors, the process of weight updating serves to train a network. If a 
peu-adigm is capable of using the input vectors’ features to converge to a network 
transfer function that accurately distinguishes each class, the network has success¬ 
fully learned the problem. 

To assess the degree of learning achieved by the networks on each problem, a 
test loop sequentially accessed all designated test vectors from the data set. The test 
vectors were originally segregated from the training vectors and were thus “unseen” 
by the network during training. This was required to truly determine the generalizing 
capability of the trained network. 

During testing, each test vector’s components and higher-order combinations 
were fed once through the network algorithm. This entailed multiplying each input 
by its corresponding weight factor stored in memory, summing the results, adding 
corresponding neuron threshold values, and feeding the summation through the sig- 
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mold function of each artificiaJ neuron. This process was mathematically identical to 
determining the output, y, for each neuron from the exemplar vectors during train¬ 
ing. However, in this case, there was no subsequent updating of weights to further 
train the net. During testing, the resulting outputs of the network for all test vectors 
were tallied to determine the number of correct class selections, yielding a percent 
accuracy for the present state of the network (determined by its stored weights). 
Following a testing session, which culminated with a stored value for network accu¬ 
racy, the program returned to the training loop to continue updating the weights 
based on additional exposures to exemplar vectors. 

For testing purposes, the criteria for determining a correct network response to 
a test vector required that the in-class node output a value greater than 0.8 and all 
other nodes (corresponding to different classes) output less than 0.2. This criteria 
is based on the networks’ estimation of the conditional probability densities (25). 
Note that for the employed sigmoid function the output range was between 0 and 1. 
Network accuracy was calculated by dividing the total number of test vectors into 
the number of correct network responses. 

Each full network run yielded a series of accuracy values for training iteration 
intervals, reflecting network learning as a function of exposures to training data. 
The software programs averaged the results of ten separate learning runs of each 
network on each problem to create the plotted results shown in Chapter IV. For 
the third-order network operating on a statistically normalized version of Ruck’s 
moment data set, results were averaged from only five runs to limit the extensive 
computation. Averaging the results from several runs served to smooth any unique 
convergence results caused by a specific combination of initial random weights for 
any one training session. 
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S.2 Development of an Image Classification Network 

Several versions of a biologic 2 Jly-motivated image classification network were 
developed and tested in the second phase of this thesis. The network algorithm was 
inspired by the biological research observations of cortical fimctioning detailed in 
Chapter II and, in part, by research conducted by Le Cun et al (16). The network 
operated on a data set consisting of segmented pixel arrays of forward-looking in¬ 
frared (FLIR) images of tanks, trucks, target boards, and clutter. The first layer of 
the network essentially correlated each image with four Gabor functions of differing 
angular orientations. This process was intended to emulate, in a highly simplified 
model, the stimulus-specific responsiveness of biological neurons in the visual cortex 
to varying spatial frequencies, or textures. Chapter II discussed neurophysiological 
evidence of cortical neurons functioning ais Gabor function detectors. 

The second layer of the image classification network performed a local averaging 
routine on the outputs of each Gabor “orientation column.” In a simplistic way, this 
process mathematically modeled phase synchronous firing of neurons in the visual 
cortex, an ol -served biological phenomena discussed in Chapter II. This higher layer 
was employed to gleeui larger global properties from am input image, rather than 
smaller stimulus details and noise extracted in the first layer of the network. 

Processed outputs from the second layer, serving ais extracted features of the 
input pixel data, were fed to the final output layer of the network for classification. 
The output layer consisted of four sigmoidal nodes, each representing one of the 
four desired output classes. Subsequent versions of the network modified the final 
output layer to find a preferred claissification scheme. Two network variations em¬ 
ployed high-order networks at the output layer to determine if various second-order 
combinations of extracted features could successfully separate the data. Another 
version implemented a multilayer perceptron with a vau-iable number of hidden layer 
nodes. Final network versions statistically normalized the features extracted from 
the Gabor wavelet correlations and local averaging routine, feeding the results to 
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the output classification layer. A diagram of the basic feature extraction algorithm 
(without details of each output layer variation for each network version) is depicted 
in Figure 18. 

S.2.1 Production of the Image Data Set A data set was developed for testing 
the image classification network. The data set consisted of integer gray-scale values 
from pixel arrays representing 88 separate FLIR images. Using software developed 
and modified by Ayer (1) to window and excise pixel arrays of images contained 
within larger scenes, 21 tanks, 23 trucks, 23 tcirget boards, and 21 clutter images 
were segmented from FLIR scenes used by Roggemann in earlier work at AFIT. 
Ayer’s software was run on VMS machines in the AFIT Graphics Laboratory to 
extract the 88 images zind convert each to a list sequence of integer gray-scale pixel 
values ranging from 0 to 255. Each image was windowed within a rectangular pixel 
array of 63 rows by 128 columns, yielding 8064 integer component values per image. 
Images were selected to minimize variations in size and were roughly centered and 
positioned upright within the extraction window. 

Each of the 88 segmented images was assigned a unique index number, from 0 
through 87, inserted before the first of the 8064 pixel veilues comprising each image’s 
list, A class identification number was then appended to the taiil of each image’s 
numerical list to identify the image class: 0 for target board, 1 for clutter, 2 for tank, 
or 3 for truck. Each image’s completed sequence of numbers consisted of 8066 integers 
including an index number, 8064 gray-scale pixel values, and a class identification 
number, in that order. The 88 separate numerical sequences representing the 88 
segmented images were finally ported to one of AFIT’s ELXSI 6400 computers and 
concatenated into one large file comprising the data set. The numerical d la set for 
the FLIR images was thus structured in a “vector-style” format for insertion into 
the network software. 
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Figure 18. Biologically-Motivated Image Classification Network Algorithm 
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S.2.S The Biologically-Motivated Image Classification Algorithm The soft¬ 
ware code for the image classification algorithm was written in C programming 
language aad compiled and run on two ELXSI 6400 computers at AFIT. Programs 
were written allowing the user to readily vary the number of exemplar versus train¬ 
ing images designated from the input data set, the number of training iterations 
rim between network accuracy tests, the total number of training iterations, and 
the learning rate, 77 . Based on the FLIR image data set produced for this thesis, 
the number of numerical components per input image was hard set at 8064. The 
number of exemplar vectors and test vectors were both set at 44. Lezurning rates 
used for various program runs were either 0.35 or 1. The software did not store and 
average results from multiple runs due to the computationally intense nature of the 
algorithm and time constraints on the project. 

Program arrays for the index, class, and image gray-scale pixel values were first 
loaded from the data set. For each program run, the random function was seeded 
from the com.puter’s time clock. Arrays for weights and thresholds were filled with 
random floating point values between -0.5 and 0.5, as values in this range fall within 
the region of greatest change for the sigmoid function. To preclude saturation of 
the sigmoid function, these weights and thresholds, along with the learning rate rj, 
were subsequently divided by the node fan in, counteracting the “explosive” eft'ect of 
huge numbers of node inputs. Also, as detailed in a separate section of this chapter, 
four arrays of constrained hard-wired weights for the first processing layer were filled 
with 2-D discrete Gabor functions. 

Training and testing loops contained identical network propagation routines. 
The training loop, randomly selecting image exemplars for network propagation, 
additionally contained standard backprop»agation routines for updating weights. The 
testing loop sequentially accessed designated test images, all “unseen” by the network 
during training, and determined if the network’s propagated outputs matched the 
desired classes. Correct test outputs were tallied and divided by the number of 
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test images, yielding a percentage for network accuracy. A correct network output 
Tvas defined by an in-class node output greater than 0.8 and all other node outputs 
less than 0.2 (for the output layer). Printouts for eeich network run displayed a 
series indicating the cumulative number of training iterations and the corresponding 
network accuracies. Specific details regarding the software are contained in Appendix 
B. 

The proposed propagation algorithm was intended to emulate, in a simple way, 
observed cortical phenomena that might relate tc the processing of visual information 
in the brain. Each input image was “blocked ofP into 8 by-8 pixel arrays to mimic 
the limited receptive fields processed by neurons in the visual cortex. Using a 50% 
overlap in both image dimensions among these “receptive fields,” an input image 
of 63 rows by 128 columns of pixels was segregated into 434 overlapping “receptive 
fields” (14 rows by 31 columns). Each 8-by-8 pixel “receptive field” was processed 
by a separate node in the first processing layer of the network. This was meant to 
imitate a recepti ve field “viewed” by a biological neuron within a cortical orientation 
column of the /isual cortex. 

An engineering judgement choosing the 8-by-8 size of the segregated pixel 
blocks was made based on the observed size of physical features in images from the 
data set. Ideally, the “receptive field” sizes were intended to capture local texture 
changes within the different images that might aid the network in discerning the 
classes. Also, due to the computationally intense nature of the proposed network, the 
50% overlap among “receptive fields” was chosen to limit the number of processing 
nodes employed at higher levels of the adgorithm. This 50% overlap could be used 
as the “mother wavelet ’ for a multiresolution hierarchy. 

Processing an input image through the first processing layer of the network 
entauled the correlation of each separate “receptive field” with four discrete Gabor 
function wavelets. A correlation multiplied each of 64 constrained hard-wired weights 
representing a Gabor wavelet by the corresponding pixel values in a ‘^ receptive field” 
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and summed the resulting 64 products. Details regarding these Gabor function 
representations of the constrained weight sets are discussed in a later section of 
this chapter. It is also noteworthy that each “receptive held” was first Lambertized 
and contrast normalized prior to correlation with the Gabor functions. Details of 
the contrast normalization routine, along with justification for employing it in the 
algorithm, sure also covered in a later section of this chapter. 

The correlated result of each “receptive field” with each Gabor function was 
fed to a unique sigmoidal node. The network structure of this processing layer thus 
consisted of 434 nodes for each of the four Gabor functions. Grouping 434 nodes, one 
per “receptive field”, for each Gabor function mimics the organization by “orientation 
columns” of neurons in the visual cortex. From this layer, the sigmoid outputs from 
each Gabor orientation group fed a subsequent layer of nodes, also segregated into 
four groups to separately process each Gabor orientation. This subsequent layer of 
nodes, linear rather them sigmoidal functions, served to emulate the phenomena of 

1- .^ase synchronization of neuronal firing discussed in Chapter II. 

The phase synchronization layer, with four groups of nodes uniquely process¬ 
ing the corresponding four Gabor orientation groups, contained 270 nodes per group. 
Each node in this layer served to average the outputs of 5-by-5 blocks of nodes from 
the Gabor function layer below. Just as the Gabor correlation layer preserved the 

2- D locational nature of the image input data in its node arrays, the phase syn¬ 
chronization layer also preserved the 2-D locationad information in processing 5-by-5 
node blocks from the Gabor correlation layer. The 5-Dy-5 node blocks from the 
Gabor correlation layer, with meocimum overlap (shifting one node in each dimen¬ 
sion), yielded 10 rows by 27 columns of blocks per Gabor group. Each linear node 
in the phase synchronization layer averaged a separate block of 25 nodes by comput¬ 
ing the arithmetic mean of their outputs. The resulting layer structure consisted of 
four groups of 270 nodes per group, each node computing a unique average output 
from the layer below. This created 1080 output values representing extracted fea- 
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tures transformtd from the original pixel data for each image. These 1080 numerical 
outputs (per Image) served as inputs to the final classification (output) layer of the 
network. Through its local averaging routine, the phase synchronization layer was 
intended to “bind” local details into larger, global features of an image. Figure 18 
illustrates the propagation hierarchy through the network layers. 

The final output layer of the network consisted of four sigmoidal classification 
nodes representing target board, clutter, tank, and truck. For each of these “top” 
nodes, outputs greater than 0.8 were considered high (in-class) emd outputs less 
than 0.2 were considered low (out-of-class). Several different versions of the output 
classification layer were tested. These network variations were needed to fully assess 
the capability of the biologically-motivated feature extraction algorithm. 

3.2.2.1 Variations of the Output Classification Layer The first network 
version employed a single layer of nodes for output classification. The output layer 
consisted of four sigmoidal nodes fully connected to nodes from the phase synchro¬ 
nization layer. This configuration assumed linear separability of the features ex¬ 
tracted from the images’ pixel data. 

In following the first part of this thesis work investigating high-order networks, 
two subsequent versions of the image classification network implemented high-order 
classification schemes. One of these versions multiplied second-order combinations of 
the 270 node outputs within each of the four separate orientation groups of the phase 
synchronization layer. The resulting 145,340 total multiplied combinations were fed, 
along with the unmultiplied features, into all four output classification nodes in 
the single-layer perceptron. The other high-order scheme multiplied second-order 
combinations of corresponding node outputs, or locations, among the four groups cf 
the phase synchronization layer. With 2700 multiplied second-order combinations, 
this approach proved far less computationally intense than the previous high-order 
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network version. The high-order network versions assumed separability of second- 
order combinations of extracted feature values. 

For more thorough testing, the network was modified to employ a multilayer 
perceptron, with two layers of traiined weights, as the output classification layer. The 
hidden layer of the multilayer perceptron allowed foi a variable number of nodes while 
the ultimate output layer consisted of four nodes, one representing each class. Vary¬ 
ing the number of hidden layer nodes provided more testing versatility by altering 
the order of mathematical computation of the classifier. The intent was to variably 
scale the computational complexity of the network to eventually achieve solution 
convergence. Software code implementing the various network versions is included 
in Appendix B. 

Noting the initial failure of these classification schemes in achieving solution 
convergence on the extracted features, it was postulated that these features might 
require statistical normalization to enable successful classification. Following the 
Gabor wavelet correlations and the local averaging algorithm, a vertical normaliza¬ 
tion routine was used to normalize, by column, the extracted feature components, 
X,-, of all 88 image vectors. The mean, p, and standard deviation, cr, for eacl fea¬ 
ture component column were calculated using only the exemplar vectors. All image 
vector components were then vertically normalized using the equation (21): 

X {{normalized) —- (13) 

<7 

Feeding the statistically normalized feature values into the multilayer perceptron, the 
most versatile of the classification schemes, the network exhibited learning with in¬ 
creased training iterations. The statistically normalized features were also fed to the 
high-order network version multiplying second-order combinations of corresponding 
node positions among the four orientation groups. This high-order network ver¬ 
sion also exhibited learning on the normalized features. Detailed result of network 
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testing are covered in Chapter IV. 


3.2.2.2 Lambertization and Contrast Normalization of Pixel Values To 
ensure network learning of relative differences among pixel values, rather than of 
absolute pixel values in an image, a normaJization routine was implemented to op¬ 
erate on 8-by-8 grouped blocks of pixel values prior to network processing. The 
8-by-8 pixel array blocks were those sectioned off as “receptive fields” for the Gabor 
function detectors in the first network layer. To simplistically emulate biological 
retinal preprocessing, normalization of the gray-scale pixel data was accomplished 
using Lambertization and contrast normalization algorithms (13, 14). 

Employing each 8-by-8 “receptive field” as a 2-D windowed region of pix¬ 
els within the image, Lambertization was accomplished by computing the average 
brightness. A, of the C4 pixel values within each window such that: 

I 64 

A=-.i:x, (14) 

where x,- represents each pixel value in the 8-by-8 window. A local contrast, xf^, 
was then computed for each of the 64 pixel values by finding the difference between 
the average window brightness and each pixel value: 

xf^ = Xi — A (15) 


The Lambertization algorithm was used to preserve local changes in image intensity 
while eliminating systematic brightness v^lriations. (13, 14) 

Following Lambertization, contrast normalization was used to mathematically 
convert vector length (“energy”), E, to unity. First, vector length Wcis calculated by: 


E = 


64 



(16) 
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Each normalized pixel value, z{{normalized), was then calculated by: 


Xi{normalized) 



(! 7 ) 


This contrast normalization procedure was used to enh 2 mce contrast where it was 
too low within the image and reduce contrast where it was too high. (13, 14) 


S.2.2.3 Discrete 2-D Gabor Functions Comprising First Layer Weights 
The network’s first layer of nodes, consisting of four gi’oups of 434 sigmoids per 
group, processed correlations of four unique Gabor functions with each of the 434 
“receptive fields” sectioned off within each image. The 2-D Gabor functions chosen 
for implementing the first processing layer of the network were 8-by-8, discretized 
wavelets possessing the same periodicity but characterized by different angular ori¬ 
entations (0, 45, 90, and 135 degrees). The 64 values corresponding to each discrete 
Gabor function served as constrained, hard-wired weight factors for inputs feeding 
the sigmoids within each Gabor orientation group. The four discrete Gabor functions 
were mathematically generated in software to fill four arrays of constrained weights 
with values between -0.5 and 0.5. This limited the Gabor function weight values 
to the region of greatest change for the sigmoid function. The following well-known 
Gabor equation (14), with selected cr values, was used to generate each function, 
9{x,yJx,fv)- 

9{x,yjxjy) = 0.5 *exp[-7r (y + y)] * cos[2 ♦ tt * (/^ y)] (18) 

where parameters fx and fy represent spatial frequencies covering two dimensions 
and input variables x and y were each discretely varied from -0.7 to 0.7 in increments 
of 0.2 to fill the 8-by-8 Gabor function weight arrays. Figure 19 depicts a plotted 
example of a Gabor function generated in this manner. 

The four different Gabor orientations were uniquely produced by varying the 
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Figure 19. Plotted Exain^de of a Mathematically Generated 8-by-8 Discrete Gabor 
Function 

combination of the fx and parameters inserted in the function. For a wavelet of 
0 degree orientation, /* = 1 and fy = 0, whereas a wavelet of 90 degree orientation 
was produced by /* = 0 and fy ~ Wavelets at 45 degrees and 135 degrees, 
possessing the same spatial periodicity as the 0 degree and 90 degree wavelets, were 
generated by combinations of = 1.4142 with fy = 1.4142, and fx -■ 1.4142 with 
fy = —1.4142, respectively. 

This chapter covered the details of the proposed high-order and image classi¬ 
fication neiwork algorithms. The following chapter discusses the testing results and 
performance of these networks. 







IV. Results 


4.1 High-order Network Performance Results 

Chapter III covered details of the algorithms for the second and third-order 
network classifiers. The following plots depict the learning performance results for 
both of these networks applied to the three test classification problems. Results for 
the second and third-order networks, each tested using two different tj values, were 
plotted separately. Baised on previous neural network research at AFIT, tj = 0.35 and 
rj = 1.0 were selected as learning rates for testing each network (21). The following 
plots show the averaged results from multiple runs for each network on each test 
problem. 

4.1.1 2-D XOR Problem Network learning performance results for the XOR 
data are depicced in the following plots. While both networks eventually converged 
to 100 % accuracy, the second-order network outperformed the third-ordei network, 
reaching maximum accuracy in far fewer training iterations for a given 7 . Also, for 
either network, a greater g yielded moie rapid solution convergence. 

It is noteworthy that, for the XOR problem, the higher-order networks did not 
seem to outperform standard feedforward, multilayer perceptrons trained through 
backpropagation. Research by T^lrr shewed that, for the XOR problem, a multi¬ 
layer perceptron network achieved 100% accuracy in 1300 training iterations. Tarr’s 
algorithms employed sigmoidal nodes and a learning rate of 7 / = 0.3 for backpropaga¬ 
tion training. Tarr’s criteria for correct network class selections required an output 
greater than 0.8 for the in-class output node and outputs less than 0.2 for cut of-class 
output nodes. (27) 

4.1.2 2 -D Mesh Problem The following plots illustrate network learning for 
the mesh problem. Neither network proved capable of learning the mesh test data to 
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Second“order Network Results for XOR 
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Figure 20. Second-order Network Learning XOR Data (Eta=0.35) 
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Second-order Network Results for XOR 

Eta-1; Fan In Hard Set to 1 



Figure 21. Second-order Network Learning XOR Data (Eta=1.0) 
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Third-order Network Results for XOR 
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Figure 22. Third-order Network Learning XOR Data (Eta=0.35) 
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Third-order Network Results for XOR 

Eta-1; Fan In Hard Set to 1 



Training iteration 


Figure 23. Third-order Network Learning XOR Data (Eta=1.0) 
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100% accuracy. Even after one million training iterations, the second-order network, 
for rj = 0.35 and r) = 1.0, achieved maximum accuracies of 64% and 66%, respec¬ 
tively. For both learning rates, the network actually approached approximate final 
accuracies within 500,000 training iterations. 

The third-order network yielded ne<irly identical results. Through one million 
training iterations, the network maintained 65% accuracy for rj = 0.35 and 67% 
accuracy for r/ = 1.0. As with the second-order network, the third-order network 
also converged to approximate final accuracies long before the millionth iteration 
was reached. Research by Ruck showed that a multilayer perceptron network trained 
via backpropagation could achieve nearly 70% accuracy on the mesh data in 8000 
training iterations, outperforming the higher-order networks investigated (24:61). 
The following plots convey the convergence properties of the higher-order networks 
and their limited maximum learning accuracies for the mesh problem. 

4 . 1.3 22-D Moments of Pixel Data (Ruck Data) The following plots depict 
network learning performance on Ruck’s twenty-two-dimensional computed moments 
from pixel data of military vehicle images. As with the mesh data, this problem 
proved beyond the learning capabilities of both high-order networks, neither achiev¬ 
ing 100% accuracy. The second-order network, for learning rates of 0.35 and 1.0, 
maintained approximate maximum accuracies of 31% and 33%, respectively, through 
one million training iterations. The third-order network maintained these same fi¬ 
nal accuracies through 500,000 training iterations. As depicted in the result plots, 
both networks actually converged to their maximum learned accuracies within a 
few hundred-thousand training iterations. Tarr’s research showed that a multilayer 
percentron network (with 77 = 0.3) eichieved approximately 75% accuracy in 50,000 
training iterations on Ruck’s data (27). It wets apparent that the higb-order mul¬ 
tiplied combinations of the moment features were not suitable for separating the 
classes. 
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Second-order Network Results for Mesh 
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Figure 24. Second-order Network Learning Mesh Data (Eta—0.35) 
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Second-order Network Results for Mesh 
Eta-1; Fan In Hard Set to 1 
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Figure 25. Second-order Network Learning Mesh Data (Eta—1.0) 
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Third-order Network Results for Mesh 
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Fig>are 26. Third-order Network Learning Mesh Data (Eta=0.35) 









Third-order Netv/ork Results for Mesh 
Eta*1; Fan In Hard Set to 1 



Figure 27. Third-order Network Learning Mesh Data (Eta=1.0) 
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Second-order Network Results for Moments 
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Figure 28. Second-order Network Learning Ruck Data (Eta=0.35) 
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Figure 29. Second-order Network Learning Ruck Data (Eta=1.0) 
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Third-order Network Results for Moments 
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Figure 30. Third-order Network Learning Ruck Data (Eta=0.35) 
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Figure 31. Third-order Network Learning Ruck Data (Eta=1.0) 
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4.S Imare Classification Network Performance Results 

Chapter III discussed details of the proposed image classification algorithm and 
the implemented variations of its output classification layer. Testing results for the 
initial version of the network, employing a single layer of percepirons as the output 
classification layer, indicated that this scheme was unable to separate the features 
extracted (without statistical normalization) from the Gabcr correlations and phase 
synchronous summations. Succe.ssful learning by this network version would have 
required linearly separable extracted features. All network versions v^ere tested with 
learning rates of 0.35 and 1. 

Both high-order versions of the image classifier also proved unable to achieve 
solution convergence on unnormalized features through one million training itera¬ 
tions. As with the previous network variation, these versions employed a single 
layer of perceptrons for the output classification layer. However, the high-order net¬ 
work implementations fed second-order multiplicative combinations of the extracted 
features to the output sigmoids. Learning success for these network versions conse¬ 
quently required separability of the second-order combinations of the unnormalized 
features. 

Another network version emplo^'ed a multilayer perception with a variable 
number of hidden layer nodes at the output classification level. Network testing 
with this versatile classifier was intended to exhaust research troubleshooting pos¬ 
sibilities at the highest processing level of the algorithm. Solution convergence was 
not achieved through 200,000 training iterations, employing ten, thirty, and fifty 
nodes in the hidden layer of the perception, and using 77 = 1.0. Following the initial 
faulure of this version of the network, it w^ls postulated that extracted features might 
require statistical normalization prior to processing by the classification scheme. A 
subsequent version of this network, performing statistical normalization on extracted 
features before feeding them to the multilayer perceptron, exhibited definite learn- 
'ng with increased training iterations. The number of nodes in the hidden layer of 
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the perceptron classifier were varied using two, three, five, ten, thirty, and fifty hid¬ 
den nodes. Learning performance results from some of the single program runs axe 
depicted in Figure 32, Figure 33, Figure 34, and Figure 35. 

The statistically normalized features were also fed to a high-order network ver¬ 
sion multiplying second-order combinations of corresponding node locations among 
the four orientation groups. High-order network runs learned and converged, within 
a few thousand training iterations, to maximum accuracies less than 10%. This 
limited performance indicated that second-order combinations of the features were 
poorly suited for separating the problem classes. 

The multilayer perceptron version of the network achieved a 43.2% maximum 
accuracy, better than the higii-order classifier, yet still limited in learning capability. 
It is likely that the Gabor functions correlated with the input images did not possess 
the optimum paxameters for extracting features to distinguish the problem classes. 
An ideal set of Gabor functions, with optimum orientations, periodicities, dilations, 
and window sizes, may yield a significantly higher majcimum learned accuracy for 
the image data set. This is discussed further along with other recommendations in 
the final chapter. 

This chapter covered performance results of the neural network algorithms 
developed for this thesis project. The following chapter discusses conclusions and 
recommendations. 
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Figure 32. Image Classification Network Results: 2 Hidden Layer Nodes 
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Image Classification Network Results 
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Figure 33. Image Classification Network Results; 10 Hidden Layer Nodes 
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Figure 34. Image Classification Network Results: 30 Hidden Layer Nodes 
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Figure 35. Image Classification Network Results: 50 Hidden Layer Nodes 
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V. Conclusions and Recommendations 


The previous chapter detailed testing and performance results for the neural 
network ailgorithms developed in this thesis project- This chapter offers conclusions 
and recommendations for continued rese 2 irch in this area. 

5.1 Investigation of the High-Order Neural Network Classifiers 

The software-based algorithms for the second and third-order neural networks 
successfully implemented and tested both high-order classification pciradigms. For 
both networks, learning accuracy increased with sequential training iterations un¬ 
til reaching a maximum accuracy for each problem. For the three classification 
problems tested in this research, the second-order network generally converged to 
a maximum learned accuracy in fewer trciining iterations, while both high-order 
networks achieved the same approximate maximum learned accuracies for the test 
problems. Consequently, for the classification problems investigated, the third-order 
network seemed to offer no advantage over the second-order network. Furthermore, 
the third-order network was more computationally intense and required significantly 
more computer run time than the second-order network. However, further research 
may prove the third-order network more capable for solving some ciaasification prob¬ 
lems. 

Both networks converged to a learned accuracy of 100% for the XOR. problem, 
the second-order network learning the pioblem much faster (in fewer iterations) 
than the third-order network. The high-order networks achieved maximum learned 
rxcuracies of approximately 66% on the mesh problem and approximately 32% on the 
Ruck data. Increasing the learning rate, g, merely reduced the number of training 
iterations required for solution convergence. 

Developing high-order algorithms beyond the third-order may produce net- 
wuikc more capable of solving these classification problems, although each added 
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order vastly increases computational requirements. This is due to explosive growth 
in the number of cross-product ajid exponential combinations of data features that 
must be calculated and subsequently fed to the sigmoid functions in each node. Also, 
it is apparent that a high-order network, of a specific order, is limited with respect 
to the classification problems it can solve. 

A lesson learned during the early developmental stages of the high-order al¬ 
gorithms was the need to normalize all inputs fed to the artificial neurons. These 
included the higher-order multiplied combinations of the data feature values (cross 
products and expcnentials) as well as the data set features. If data set feature values 
are relatively small (between 0 and 1), their higher-order multiplied combinations 
are even smaller. Miniscule high-crder inputs may not contribute to weight updat¬ 
ing, thus precluding network learning. Prior to normalizing the high-order inputs, 
network learning accuracies widely vacillated even after many training iterations, 
never settling to a maximum learned accuracy. I hese e-ratic fluctuations in network 
accuracy were eliminated following insertion of normalization code modules for the 
high-order network inputs. 

The high-order networks did not outperform feedforward multilayer percep- 
trons developed by Ruck and Tarr and tested on the same classification problems. 
In selecting a general classification network for learning a data set, the multilayer 
perceptron may prove superior in many cases. The multilayer perceptron allows for 
.Tying the number of nodes in its hidden layers to “match” the order of a problem 
and achieve maximum learning. A good, easily modified classifier should be able to 
separate output classes for many different problems, even with no a priori knowledge 
of the organization of the input feature data. For example, examining a plot of the 
segregated class regions of the 2-D mesh data (Figure 16), it seems likely that using 
radial basis functions or logp, ^ inputs in a network might readily solve this specific 
problem (13). Yet while a priori knowledge of a specific problem may be highly 
useful in developing a network solution, such information is not always available, 
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particularly for data sets consisting of higher dimensional vectors. 

5.2 Investigation of the Image Classification Algorithm 

The inability of the initial version of the image classification network to achieve 
solution convergence prompted the development of subsequent variations in the out¬ 
put classification layer. These network variations, which included two high-order im¬ 
plementations and a multilayer perceptron with a variable number of hidden nodes, 
were employed to thoroughly test and verify inadequacy of the feature extraction 
scheme by exhausting possibilities of deficiencies in the output classification layer(s). 
After employing statistical normalization on extracted features prior to classification, 
the multilayer perceptron and a high-order network were capable of converging to 
solutions as learning followed increased training iterations. 

To extract data features more suitable for distinguishing designated output 
classes, future optimization of the algorithm may vary aspects of the feature ex¬ 
traction scheme, including the discrete wavelet correlations and phase synchronizing 
local averaging routine. Worthy of further exploration would be the use of discrete 
Gabor functions characterized by different periodicities, orientations, dilations, and 
window sizes relative to the four Gabor wavelets employed in this research. Vary¬ 
ing these aspects in the correlation process would extract different features from 
the image pixel data. Tarr proposed a methodology in his draft dissertation for se¬ 
lecting optimum Gabor function parauneters for image feature extraction (27). The 
phase synchronization scheme may also be modified to implement inhibition among 
the four phase synchronized groups derived from the four different Gabor function 
correlations. Successful future algorithms could implement an energy thresholding 
scheme in the Lambertization and contrast normalization routine, possibly optimiz¬ 
ing the program (13). Varying the degree of overlap for pixel windows in the input 
plane and node blocks processed by the local averaging routine may also serve to 
optimize the feature extraction algorithm. The classifier may ultimately be modified 
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to lecognize problem classes with position, scale, and rotation invariance. 

This ie.search project employed ELXSI 6400 computers, with single network 
programs often running for days, and sometimes weeks. Due to the computational 
intensity of these image classification algorithms operating on pixel data, it may be 
useful in the future to run software models on a Cray computer. Furthermore, care 
should be teiken to avoid exceeding machine memory capacity when running intensive 
processes. This problem is exacerbated when using double precision variables to 
preserve critical small differences in interim program values. 

A recommendation that may improve the tcdlying of network convergence en¬ 
tails the use of less stringent criteria for defining correct network class responses 
during testing. This thesis research defined correct network responses from the out¬ 
put layer as the in-cl 2 u.s node firing greater than 0.8 and all other nodes (out-of-class) 
firing less than. 0.2. A less exacting requirement for determining a correct class re¬ 
sponse might be simply declaring the highest output node as the in-class node. 

The FLIR image data set, consisting of tanks, trucks, target boards, and clut¬ 
ter, proved to be a challenging classification problem. Reducing the output classes, 
or running the image classification network on a more easily separable data set, maj'^ 
also facilitate optimization of the algorithm. 
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Appendix A. High-order Network Code 


78 


/* single-layer, second-order neural network algorithm with 
automatic vertical normalization of all input combos. 
Results averaged over 10 runs. */ 

/* Name; N0RMNET2.C */ 


#lnclude <stdio.h> 
#lnclude 'diath.h> 


#deflne numneur 4 

#deflne numexvec 58 
#define ntimttvec 23 
#deflne numfeat 22 

#deflne niunzuns 1000 
#deflne numtm 1000 


/* Specify total number of neurons across slab */ 
/* Number of neurons « number of classes */ 

/* Specify number of exemplar vectors */ 

/* Specify number of test vectors */ 

/* Specify number of features per vector */ 

/* numfeat also represents # FEATURE COLUMNS */ 

/* Specify number of result plot points */ 

/* Specify # training Iters between test/plot */ 


^define eta 1.0 


/* Specify training factor */ 


char *malloc(); 
maln(){ 


/* Must define for use of malloc thru program */ 
/*Need for this statement is compiler-dependent*/ 


/* flag[] holds ID number of each data vector, class[] holds training 
class associated with each data vector, x[][] holds feature 
components of data vectors, quadln[][] holds second-order multiplied 
combos of X values, wlin[](] and wquad[](] hold weights for all inputs 
(including higher order combos), theta[] holds threshold values for 
each, neuron */ 


int *flag; 

Int 1, j, n, randnum, n\imquad, qcount, fanln; 

Int row, col, avgrun; 

int neurtcnt, numrlght, blgloop, tmloop, testvec, Itnum; 
float d, y, error, neweta, accuracy, addup; 
float *x, *class, *wquad, *quadln; 

float theta[niunneur], wlin[numneur*numfeat], avrun[10][numruns]; 
float flnalavg[n\imruns] ; 

float siun, llnsum, quadsum, mean, s, sumsqdlf; 

/* This module opens an input file and reads data into the 
arrays flag[], class[], and x[][] */ 

FILE *fp; 

fp“fopen("ruckdata","r"); /* Insert correct data file name here */ 

if(fp—NULL) exit(O); 

/* D3mamlc alloc of arrays prior to loading in data file values */ 

flag - (int *) malloc((numexvec+numttvec) * sizeof(int)); 
class - (float *) malloc ((niimexvec-fnumttvec) * sizeof (float)) ; 

X (float *) malloc((numexvec+numttvec)*numfcat * sizeof(float)) ; 
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for(i-0; KCnumexvec+numttvec); ++i)( 
fscanf(fp,"%d",&flag[i]); 

for(j-0; J<numfeat; ++j){ 
f scanf (fp, "%f " ,&x[ i*nuinf eat+j ] ) ; 

} /* end inner loop */ 

fscanf(fpc "%f",&class[1]>; 

) /* end outer loop */ 

fclose(£p); 

/* This test module prints out input data containing first two and 
final feature components of each vector; checks the above file 
opening module (comment out for final program) */ 

/* for(i-0; i<(ntimexvec+numttvec) ; ++i) { 

printfC'For flag-%d, class-%f, check features are %f %f %f\n", 
flag[l], class[i], x[i*numfeat+0], x[i*numfeat+l], 
x[i*numfeat+(numfeat-l)]); 

> */ 

/* This module establishes numbers of second-order 

combos of feature inputs FOR EACH VECTOR (without commutative 
redundancy). Will be used for array allocation. */ 

n\imquad-0; 

for(i-«l: i<-n\imf eat: i++) { 
numquad>ni\imquad+l; 

) 

/* Fill 2-D arrays for second order input combos; 
quadin[numexvec+n\imttvec] [numquad] */ 

quadlno(float *) malloc((numexvec+numttvec)*numquad*siteof(float)) 
for(row-0; row<(numexvec+numttvec); row++){ 
qcount-0; 

for(i-0; Knumfeat; i++)( 
for(j-i; j<numfeat; j++){ 

quadin[row*n\imquad+qcount]-xtrow*numfeat+i]*x[row*numfeat+j]; 
qcount-qcount+1; 

) 

) 

) 

/* Test printout module */ 

/* for(row-0; vow<(numexvec+numttvec); row++){ 
for (col-0; coKnumquad; col++) { 

printf("For row»%d col-%d quadin-%f\n",row,col, 
quadin[row*nufflquad-fcol]) ; 

) 

) */ 
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/* Begin VERTICAL NORMALIZATION of x, quadin arrays */ 

/* Vertically normalize x[n\imexvec+numttvec][numfeat] */ 
for'col"0; col<numfeat; col++){ 
sum-0.0; 

for(row-0; rov<numexvec; row++){ 
sum-sum+x[row*niimfeat+col]; 

) 

mean-sum/(float)(numexvec); 
sumsqdif-0,0; 

for(row»0; row<numexvec; row++){ 

sumsqdif-sumsqdif+(x[row*numfeat+col]-mean)*(x[ro%r*numfeat+col]-mean) 

) 

s-(float)(sqrt(sumsqdif/(float)(numexvec))); 
for(-:ow-0; row<(numexvec+numttvec) ; row++)( 
x[row*nvimfeat+col]-(x[row*numfeat+col]-mean)/s; 

) 

) /* End column incrementing loop */ 

/* Test printout loop */ 

/* for(row-0; row<(numex’.' c+numttvec); row++)( 
for (col-0; coKnumfeat; colx+){ 

printfC'For row-%d col-%d nomxvalue-%f\n" ,row,coI, 
x|.row*n\amfeat+col]); 

) 

) */ 

/* Vertically normalize quadin[numexvec+numttvcc)[numquad] */ 
for (col-0; coKnumquad, col++){ 
sura-0.0; 

for(row-0; row<numexvec; row++){ 
sum-sum+quadin[row*numquad+col]; 

) 

mean-sum/(float)(numexvec); 
s\imsqdif-0.0; 

for(row-0; row<numexvec; row++){ 
sumsqdif-sumsqdif-f 

(quadin[row*numquad+col]-nean)*(quadln[row*numq\iad+col]-mean); 

> 

s-(float) (sqrt(sumsqdii;/(float) (numexvec))); 
for(row-0; row<(numexvec+numttvec); row++){ 
q\iadln[row*numquad+coJ]-(quadin[row*numquad+col]-mean)/s; 

} 

) /* End column incrementing loop */ 

/* Test printout loop */ 

/* for(row-0; row<(numexvec+numttvec); row++){ 
for (col-0; col-Oiumquad; col++){ 

prlntf("For row-%d col-%d normquadvalue-%f\n",row,col. 
quadin [rov«rnumq\iad-«-col} ) ; 

) 

) */ 

/* End of vertical normalization */ 
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£or(avgrun-0; avgrun<10; avgrun++){ /* 10 runs thru program for avg */ 
srandom((unslgned)tlme(NULL)); /* Inlt random() seed off clock */ 

/* Loop to Initialize theta array (one for each neuron) with random 
floating point value between -0.5 and 0.5 */ 

for(l“0; Kntimneur; !++){ 

theta[i]-(float) (randomO % 101)/100.00 - 0.5; 

/* printf(’'For l-%d theta«%f\n" , 1, theta f 1]) ; */ 

) 

/* Loop to Initialize linear input weights array wlin[numneur][numfeat] 
with random floating point values between -0.5 and 0.5 */ 

for(n-0; noiuameur; n++){ 

/* printfC’For neuron # %d\n",n); «/ 
for(1-0; i<n\imfeat; !++){ 

wlln[n*n\imfeat+i]-(f loat) (randomO % 101)/100.00 - 0.5; 

/* printfC'For i-%d wlin»%f\n" . i.wlin[n*nuinfeat+i] ) ; */ 

) 

) /* end of nested loop */ 

fanin - numfeat + ntimquad; 

/* fanin - 1; */ 

neweta - (float) eta / (float) fanin; 

/* printf("neweta- %f\n",neweta); */ 

/* Loop to initialize quadratic input weights array 

wquad[numneur][numquad] with random floating point values between 
-0.5 and 0.5*/ 

wquad - (float *) malloc(n’imneur*numquad*sizeof(float)); 
for(n-0; n<n\imneur; n++){ 

/* printfC'For neuron # %d\n",n); */ 
for(i-0; Knumquad; i+*-){ 

wquad[n*numquad+i]-(float)(randomO % 10i)/100.00 - 0.5; 

/* printfC'For i-%d quadweight"%f\n",l,wquad[n*numquad+i]); */ 

) 

) /* end of nested loop */ 

/* End of initialization; Begin loops for training and testing */ 
for(bigloop-l; bigloop<-numruns; bigloop++){ 

/* Begin training loop */ 

for(tmloop-0; tmloop<n\imtm; tmloop++){ 

/* Randomly select an exemplar vector row (flag #) */ 
randnum-(randomO % niunervec) ; 

/* printf("Randomly selected vector flag number is %d\n".randnum); */ 

for(n-0; n<numneur; n-M-){ /* Loop to update weights for each neuron *'</ 

llnsum-0.0; 

for(i-0; Knumfeat; !++){ 

1 insum-1 Insum+wl 1 n [n*numf eat+i ] *x [ randnum*nujnf eat+i ] ; 

) 
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qviadsum-0.0; 

for(l«*0; Knufflquad; 1++) { 

quadsum“quadsum+wquad[n*nuinq\iadvi]*quadin[randnum*ntimquad+i] ; 

) 

sum»*l insiyn+quadsuiB; 
if ((sum+theta[n])<-20.0) (y»0.0;} 

else if ((st«s+theta[n])>20.0) (y-l.O;) 

else {y-(float)(1.0/<1.0fexp(-(suin+theta[n]))));) 

/* Criteria for training each class to fire a specific neuron */ 
if (class[randnum]—(float) (n)) (d-l.O;) 
else (d-O.O;) 

error-y*(1.0-y)*(d-y); 

/* prlntf("error-»%f\n" .error); */ 
for(i-0; i<numfeat; i++){ 

wlin [n*n\imf eat+i ] «wl in [n*numf eat+i ] +neweta*error*x i randnum'^iumf ea t + i ) 
/* printfC'For i-^ld wlin-%f\n'*, i ,wlin[n*nvimfeat+i j ) ; */ 

) 

for(i-0; i<numqua d; i++){ 

wqua d[n*numquad+i]-wquad[n*numquad+i] + 

neweta*error*quadin[randnum*numquad+i]; 

/* printfC'For i»%d quadweight*»%f\n",i,wquad[n*nuinquad+i]) ; */ 

} 

theta[n]-theta[n]+neweta*error; 

/* printfC'For neuron # %d theta+%f\n",n,theta[nj); */ 

) /* end of weight training loop for each neuron's weights */ 

) /* end of complete training loop */ 

/* Begin test loop to run thru each test vector */ 
numrlght-O; 

for(testvec-numexvec; testvec<(nuinexvec+numttvec); testvec++){ 

/* printf("testvec-%d\n",testvec); */ 

neurtcnt-0; 

for(n-0; n<numneur; n++)( 
llnsum-0.0; 

for(1-0; i<n\unfeat; i++) { 

llnsuin-llnsvun+wlln[n*nuinfeat+lj*x[testvec*nuinfeat+i] ; 

) 

quadsum-0.0; 

f or (1-0; Knumqua d; i++) { 

quadsum»q\iadsvun+wquad[n*nujnquad+i]*quadln[testvec*nuinquad+i ] ; 

) 

stun-1 insum+quadsum; 

if ((sum+thetain])<-20.0) (y-0.0;) 

else if ((suin+theta[n])>20.0) (y-l.O;) 

else (y-(float)(1.0/(1 0+exp(-(sum+theta[n]))));) 

/* printf("y-lf class[tesl:vec]-%f\n”,y,class[testvec]); */ 
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/* Test decision criterie for identifying desses */ 

if ((y>0.8) && (clets[testvec]—(,floet)(n))) {neurtcnt-neurtcnt + 1 . ) 

if ((y<0.2) && (cless[testvec j *-i’f>o«t)(n))) (neurtcnc-neurtcnt* 1 , ) 

) /* end of neurons (n) correct test loop */ 

if (neurtcnt>-«usneur) (nuaright’Mauarightel;) 

) /* end of testing loop running thru each test vector (tescvec) */ 

itnuiB-bigIoop'*nriUBtm: 

accurecy-(f lost)nu*rlght/(f lost )n\i«ttvec ; 

/♦printf("Training iteration %d Accuracy-%f\n".Icnun.accuracy) ,'^/ 
avTun(avgrun]Iblgloop-1]-accuracy; 

) /* end of bigloop */ 

) /* end of 10-run loop for averaging */ 

/* Final averaging routine */ 
for (col-0; coKnumruns ; col++)( 
addup-0.0; 

for(row«0; row<10; row++){ 
addup“addup+avrun(row][col]; 

) /* end row cdd up loop */ 

flnalavg[col]-addup/10.0, 

printf("Training Epoch #: %d Average Accuracy-%f\n", 

(col+l)*nuintm,f lnalavg[col]); 

) /* end of column loop */ 

) /* end of main */ 
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/* single-layer, third-order neural network algorithm with 
automatic vertical normalization of all input combos. 

Results averaged over 10 runs. */ 

/* Name: N0RMNET3.C */ 

#include <stdio.h> 

^include ■(3aath.h> 

#define numneur 2 /* Specify total number of neurons across slab */ 

/* Number of neurons - number of classes */ 
#define numexvec 4 /* Specify number of exemplar vectors */ 

#define numttvec 4 /* Specify n\imber of test vectors */ 

#define numfeat 2 /* Specify number of features per vector */ 

/* numfeat also represents # FEATURE COLUMNS */ 
#deflne numrtins 300 /* Specify number of result plot points */ 

#define numtm 10 /* Specify # training iters between test/plot */ 

#define eta 0.35 /* Specify training factor */ 


char *malloc(); 
main() { 


/* Must define for use of malloc thru program */ 
/*Necd for this statement is compiler-dependent*/ 


flag[] holds ID niunber of each data vector, class[] holds training 
class es<‘ociated with each data vector, x[][] holds feature 
components of data vectors, quadin[][] holds second-order multiplied 
combos of X values, cubeln[][] holds third-order multiplied combos 
of X values, vlln[][], wquad[j[]. and wcube[][] hold weights for all 
Inputs (ino'' ding higher order combos), theta[] holds threshold 
values vor each neuron */ 

int *flag; 

Int i, J, n, randnum, numquad, numcube, qcount, cubcount, fanin; 
int startcnt, upcount, row, col, avgrun; 

int neurtcnt, numright, bigloop, tmloop, testvec, itnum; 

float d, y, error, neweta, accuracy, addup; 

float *x, *class, *wquad, *wcube, *quadin, *cubein; 

float theta[ntimneur], wlin[numneur*numfeat], avrun[10][numruns]; 

float finalavg[numr\ins] ; 

float sum, linsum, quadsum, cubesum, mean, s, sumsqdif; 

/* This module opens an input file and reads data into the 
arrays flag[], classi], and x[][] */ 

FILE ★fp; 

fp«fopen("xordata","r"); /* Insert correct data file name here */ 

if (fp—NULL) exit(O): 

/* Dynamic alloc of arrays prior to loading in data file values */ 

flag » (int *) malloc((numexvec+numttvec) * sizeof(int)); 
class - (float *) malloc((numexvec+numttvec) * sizeof(float)); 

X "• (float *) malloc((nufflexvec-fn\imttvec)*numfeat * sizeof (float)) ; 
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for(i-0; i<(numexvec+nuinttvec); ++!){ 
fscanf(fp,"%d",&flag[i]); 

for(j-0; j<nTamfeat; ++j){ 
fscanf (fp, "%f" ,&x[ i*nuinfeat+j ] ) ; 

) /* end inner loop */ 

fscanf(fp, "%f",&class[l]); 

) /* end outer loop */ 

fclose(fp); 

/* This test module prints out input data containing first two and 
final feature components of each vector; checks the above file 
opening module (comment out for final program) */ 

/* for(i“0; l<(numexvec+numttvec); ++!)( 

printfC'For flag-%d, class-%f, check features are %f %f %f\n", 
flag[i], class[i], x[i*numfeat+0], x[i*numfeat+l), 
x[i'*n\imfcat+(numfeat-l) ]); 

) */ 

/* This module establishes numbers of second-order and third-order 
combos of feature inputs FOR EACH VECTOR (without commutative 
redundancy). Will be used for array allocation. */ 

ntamquad»=0; 
numcube-0; 

for(i-l; K-numfeat; i++){ 
numquad-numquad+i; 
numcube-^umcube+numquad; 

) 

/* Fill 2-D arrays for second and third order input combos: 
quadin[numexvec+numttvec][numquad] and 
cubeln[numexvec+n\imttvec] [numcube] */ 

quadin“(float *) malloc((numexvec+numttvec)*Tnimquad*sizeof(float)) 
for(row-0; row<(numexvec+n\amttvec); row++){ 
qcount-0; 

for(l-0; i<n\imfeat; i++){ 
for(j-l: J<numfeat; j++){ 

quadin[row*numquad+qcount]—x[row*numfeat+i]*x[row*numfeat+j]; 
qcount-qcotint+l; 

) 

) 

1 

/* Test printout module */ 

/* for(row“0; row<(numexvec+numttvec); row++){ 
for (col-0; coKnumquad; col++){ 

printf("For row—%d col—%d quadin-%f\n",row,col, 
quadintrow*numquad+col]); 

) 


cube In-(float *) malloc( (niunexvec:+ntimttvec)*numcube*slzeot(f loat) ) ; 

fcr(row-0; row<(nume3i:ve>:-5-nuinttvec) ; row++)( 

cubcount-0; 

startcnt-O; 

•upcount-numf eat; 

forCi'O; Icnumfeat; !++){ 

for(j-startcnt; j<nuinqua<!; J++){ 

cubeln [row*numcv\be4-cubccunt )-x [ro%»*nujnfeat+i] *quadin [row*nuinquad+j ] 
cubcount-cubcovmt'tl; 

) 

startcnt-startcnt+upcount; 
upcount-upcount-1; 

) 

) 

/* Test printout loop */ 

/* for(row"C; row<(numexvec+numttvec); row++){ 
for(col-0; coKnumcube; col++){ 

prlntf("For row-%d col-%d cube in-%f\n‘'.row, col, 
cubcin[row*nunicuba+col ]); 

) 

) */ 

/* Begin VERTICAL NORKALIZATION of x, quadi.n, cubein arrays */ 

/* Vertically normalize x[nvunexvec+numttvec][numfeat] */ 
for (col-0; coKnumfeat; col++) { 
sum-0.0; 

for(row-0; row<numexvec; row+i){ 
sum-sxun+x[row*numfeat+col]; 

) 

mean-s\im/(f loRt) (numexvec) ; 
sumsqdif-O.0; 

for(row-0; row<numexvec; row++){ 

svmisqdif-suBisqdif+(x[row*numfeat+col] -niean)*(x[row*numf eat+col ] -mean) 

) 

s»(float)(sqrt(sumsqdif/(float)(numexvec))); 
for(row-0; row<(numex7ec+numttvec); row++){ 
x[row*numfeat+col]-(x[ro%r*numfeat+col]-mean)/s; 

) 

} /* End coltimn incrementing loop */ 

/* Test printout loop */ 

/* for(row-0; row<(numexvec+numttvec); row++){ 
for (col-0; coKnumfeat; col++){ 

prlntf("For row-%d col-%d nornixvalue-%f\n",row,col, 

X[row*numfeat+col]); 

) 

) */ 
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/* Vertically novmallKe qua'lln[numoxvec+n\unttvec] [numquad] */ 

for(col-0; col<n\imquad; col+f){ 
sum"»0.0; 

for<row-0; row<numexvec; rcw++){ 
suai-sti]n+quadlii[row*nuniquad+col] ; 

) 

mean»s\im/(f loat> (numexvec); 
sumsqdlf>0.0; 

for (row-0; row<n\aiaexvec; row++){ 
sumsqdif-suinsqdif+ 

(quadin[row*numquad-fcol] -mean)*(quadin[row^mquad-(-col] -mean) 

) 

s-(float)(sqrt(sumsqdif/(float)(numexvec))); 
for (row-0; row<(niimexvec+ntunttvec) ; row++){ 
quadln[row*numquad+col]-(qvadln[row*niimquad+col]-mean)/s; 

) 

) /* End column incrementing loop */ 

/* Test printout loop */ 

/* for(row-0; row<(niimexvec+n\imttvec) ; row++) { 
for (col-0; coKntimquad; col++) { 

printf(’'For row-%d col-%d normquadvalue“%f\n",row,col, 
quadin[row*numquad+col]); 

) 

} */ 


/* Vertically normalize cubein[n\mexvec+n\unttvec] [n\uncube] */ 

for(col-0; coKnumcube ; col++) { 
stun-0.0; 

for(row-0; row<numexvec; row++)( 
sum-svun+cubein[row*numcube+col]; 

) 

mean-sum/(float)(numexvec); 
suffisqdif-0.0; 

for(row-0; rowcnumexvec; row++){ 
sumsqdif-sumsqdif-i- 

(cubein[row*numcube+col]-meaii)*(cubein[row*numcube+col)-mean) 

) 

s-(float)(sqrt(sumsqdlf/(float)(numexvec))); 
for(row-0; row<(numexvec+numttvec); ro5/++){ 
cubein[row*numcube+col]-(oubein[row*numcube+col]-mean)/s; 

} 

} /* End column incrementing loop */ 

/* Test printout loop */ 

/* for(row-0; row<(numexvec+numttvec); row-H-){ 
for(col-0; coKnumcube; col++){ 

printf("For row-%d col»%d normcubevalue-%f\n",row,col, 
cubein[row*numcube+col]); 

} 

) */ 

/* End of vertical normalization */ 





for(avgrun«0; avgrun<10; avgr\in-H-){ /* 10 mns thru program for avg */ 
srandom<(unsigned)time(NULL)); /* Init randomO seed off clock */ 

/* Loop to initialize theta array (one for each neuron) with random 
floating point values between -0.5 and 0.5 */ 

for ( i»0; Knuameur; 1++) { 

theta[i]-(float)(random() % 101)/100.00 - 0.5; 

/* printf("For i-%d theta-%f\n",1,theta[1]); */ 

) 

/* Loop to Initialize linear input weights array wlintnumneur][numfeat] 
with random floating point values between -0.5 and 0.5 */ 

for(n-0; n<n\iameur; n++)( 

/* printf("For neuron # %d\n",n); */ 
for(i«0; Knumfeat; i++){ 

wlin[n*nuinfeat'+i]-(f loat) (randomC) % 101)/100.00 - 0.5; 

/* printfC'For l»%d wlln-%f\n",i.wlin[n*numfeat+i)); */ 

) 

) /* end of nested loop */ 

/* fanin - numfeat + numquad + numcu'..-e; */ 
fanln - 1; 

neweta - (float) eta / (float) fanln; 

/* printf("neweta- %f\n",neweta); */ 

/* Loop to initialize quadratic input weights array 

wquad[numneur][ntunquad] with random floating point values between 
-0.5 and 0.5*/ 

wquad - (float *) malloc(numneur*numquad*siTeof(float)); 
for(n-0; n<numneur; n++){ 

/* prlntf("For neuron # %d\n",n); */ 
f or (1-0; Knumquad; 1++) { 

wquad[n*numquad+i]“(float)(random() % 101)/100.00 - 0.5; 

/* printf("For i-%d quadwelght-%f\n",i,wquad[n*numquad+l ]); */ 

) 

i /* end of nested loop */ 

/* 7x)op to initialize cubic input weights array wcube[numneur][numcube] 
with random floating point values between -0.5 and 0.5 */ 

wcube - (float *) malloc(numneu?:*numcube*sizeof(float)); 
for(n-0; n<numneur; n++){ 

/* printf("For neuron # %d\n",n); */ 
for(i-0; Knumcube; i++){ 

wcube[n*n\imcube-fi]-(float)(randomO % 101)/100.00 - 0.5; 

/* printf("For l-%d cubewelght-%f\n",i.wcube(n*rumcubi+lJ); */ 

) 

) /* end of nested loop */ 


5 





/* End of Initialization; Begin loops for training and testing */ 
for(blyloop«l; bigloop<“numruns; blgloop-M ){ 

/* Begin training loop */ 
for(tmloop>«0; tmloop<numtm; tmloop++)( 

/* Randomly select an exemplar vector row (flag #) */ 
randnuffl<-(random() % ntunexvec); 

/* printf("Randomly selected vector flag niimber is %d\n".randnum); */ 


for(n“0; n<numneur; n++){ /* Loop to update weights for each neuron */ 

llnsum-O.O; 

for(i-0; Knumfeat; !++)( 

linsum'"llnsum+wlln[n*numfeat+i]*x[randnum*numfeat+i]; 

) 

quadsuffl-0.0; 

for(i“0; iCnumquad; i++){ 

qua dsum-quad sum+wqua d[n*Tiumqua d+i]*quadin[randnum*numqua d+i); 

) 

cubesum-0.0; 

for(i-0; Knumcube; i++){ 

cubes\im-cubesum+wcube fn*numcube+i]*cubeln[randn\un*numcube+i ] ; 

1 

sum-1Insum+quadsum+eube sum; 
if ((sum+theta[n])<-20.0) (y-O.O;) 

else if ((sum+thetatn])>20.0) (y-l.O;) 

else {y»(float)(l,0/(1.0+exp(-(s\im+theta[n]))));) 

/* Criteria for training each class to fire a specific neuron */ 
if (class[randnum]--(float)(n)) {d-l.O;} 
else (d-O.O;) 

error-y*(1.0-y)*(d-y); 

/* printf (‘'error-%f\n" .error) ; */ 
for(i>*0; Knumfeat; !++)( 

wlin[n*numfeat+i]«wlln[n*numfeat+i]+newete*error*x[randnum*numfeat+i]; 
/* printf("For l-%d wlin-%f\n",i,wltnln*numfeat+i]); */ 

) 

for(i-0; Knumquad; !++)( 

wquad[n*numquad+l]-wquad[n*numquad+i]+ 

neweta*error*quadin[randnum^numquadti]; 

/* printf ("For l»%d q\xadwelght-%f\n",l,wquad[n*n\iraquad+l]); */ 

) 

for(l-0; Knumcube; !++){ 

wcube[n*numcube+i]-wcube[n*numcube+iJ+ 

neweta*error*cubein[randnum*numcube+l]; 

/* printf("For i-%d cubeweight-%f\n",i,wcube[n*numcube+i]); */ 

) 

theta[n]-theta[n]+neweta*error; 

/* printf("For neuron # %d theta-%f\n",n,theta[n]); */ 

) /* end of weight training loop for each neuron's weights */ 

) /* end of complete training loop */ 
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/* Begin test loop to run thru each test vector */ 
numrlght«*0; 

forCtestvec-numexvec; testvec<(numexvec+numttvec); testvec++){ 

/* prlntf<"testvec-%d\n".testvec); */ 
neurtcnt-0; 

for (n-0; n<nuimeur; n++) ( 
llns\i»-0.0; 

for (1-0; Knumfeat; !++){ 

llnsuai-llnsum+wlin[n*numfeat+l]*x[tcstvec*numfeat+l]; 

} 

quadsum-0.0; 

for (1-0; l'<^u!nquad; 1++) ( 

quadsum-quadsuffl+wquad[n^umquad-(-l]*quadln[ testvec*numquad-«-i] ; 

) 

cubestin-0.0; 

for(l-0; Knuncube; !++)( 

cubesum-cubestim+wcube[n*numcube+ll*cubein[ testvec*nviracube+i] ; 

) 

sum-1 insum-t-q\iad sum-«-cube sum; 

If ((s\im+theta[n])<-20.0) (y-O.O;) 

else If ((s\im+theta[n])>20.0) (y-l.O;) 

else {y-(float)(1.0/(1.0+exp(-(sum+theta[n]))));) 

/* prlntf("y-%f class[testvec]-%f\n",y,class[testvec]); */ 

/* Test decision criteria for Identifying classes */ 

If ((y>0.8) && (class[testvecl—(float)(n))) (neurtcnt-neurtcnt+l;) 
If ((y<0.2) && (classitestvecjl-(float)(n))) (neurtcnt-neurtcnt+l;) 

) /* end of neurons (n) correct test loop */ 

If (neurtcnt—numneur) (numright-numrlght+l;) 

) /* end of testing loop running thru each test vector (testvec) */ 

itnum-bigloop*numtm; 

8ccuracy-(float)numright/(float)numttvec; 

/*prlntf ("Training Iteration#: %d Accuracy-%f\n" , itnvuii,accuracy) ;*/ 
avrun[avgnm] tbigloop-l]-accuracy; 

) /* end of bigloop */ 

) /* end of 10-run loop for averaging */ 

/* Final averaging routine */ 
for (col-0; coKnumruns; col++){ 
addup-0.0; 

for(row-0; row<10; row++)( 
addup-addup+avrun[rov][col]; 

) /* end row add up loop */ 

fInalavg[col]-addup/10.0; 

prlntf("Training Epoch #: %d Average Accuracy-%f\n", 

(col+l)*nuiiitm,f Inalavg [col]); 

) /* end of column loop */ 

) /* end of main */ 


7 


Appendix B. Image Classification Network Code 
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/* Image classification network. This algorithm processes input 
data sets consisting of ais index number, a class number, and 
8064 components per vector (for 63 X 128 image pixel arrays). 

Input "receptive field" window.^ are set at 8 X 8 pixels . There 
are two hidden layers feeding four output nodes representing 
tanks, trucks, target boards, and clutter. Name: IMAGENET.C */ 

#lnclude <stdio.h> 

#include <math.h> 

#include <strtict.h> 

#define numeirvec 44 /* Specify number of exemplar vectors */ 

#define numttvec 44 /■*’ Specify number of test vectors */ 

#deflne nunfeat 8064 /* Specify number of features per vector */ 

#deflne numrtms 60 /* Specify niifflber of result plot points */ 

#define nurntm 1000 /* Specify # training epochs between test/plot */ 

#deflne eta 0.35 /* Specify learning-rate */ 

#defirie pi 3.1415326535d979 

/* Structure creating 4-D weight array feeding output neuron.s : 
outnode(i).block2[j].layotwgt[k] */ 

struct tagl { 

double layotwgt[270I; 

); 

struct { 

struct tagl block2[4]; 

) outnode[4]; 

char *malloc(); 

HiainO { 

int *flag, *x, *class, neurtent, itnum, gabcount; 

int i, j, k, n, bigloop, tmloop, testvec, numright, randninn, bshlft; 
int rowstpix, colstpix, rowsthl, colsthl, index, laylposn, iay2posn; 
double rfarray[64], absum, avbright, esum, energy, lay2wgt; 
double outtheta[4], sun., fx, fy, a, b; 

double laylout[4] [434] , lay?.out[4] [270], pharray[25] , y[4) ; 
double gabor[4][64] ; 
double d[4], errorout[4]; 
float accuracy; 

/* Open and read in data file to arrays fl8g[], class[], x[][] */ 

FILE *fp; 

fp-fopen("blgimage", ’’r") ; 
if(fp—NULL) exit(O); 

flag - (int *) Eialloc((ntimexvec+numttvec) * slzeof (int)) ; 
class “ (int *) malloc((numexvec+numttvec) * sizeof(int)); 

X - (int *) malloc((nuroexvec-t-numttvec)i*fnumfeat * sizeof (int)) ; 
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for(i-0; i<(nuiii®xvec+numttvec) ; 1++) { 
fscanf(fp,"%d",&flag[i]); 
forCj-O; JcSiitimfeat; j++){ 
fscanf (fp, "%d" ,&x[i*Tiutnfeat+j ] ) ; 

) 

fscanf(fp,"%d",6iclass[l]) ; 

) 

fclose(fp); 

/* Test printout module: */ 

/* for(i“0; l<(nximexvec+numttvec); i++)( 
printf('*Flag-%d\n".flag[i]) ; 

) V 

srandom((unsigned)time(NULL)); /* Init random() seed off clock */ 

/* Initialize output theta and weight arrays */ 
for(i-=0; i<4: i++) { 

outtheta[i]-(double)(randomO % 101)/100.00 - 0.5; 

) 

for(i-0; 1<4; i++){ /* Normalize initial values by fan in */ 

outtheta[i]“Outthetafi]/1080.00; 

) 

for(i-0: KA; i++){ 
for(j-0; j<4: j++){ 
for(k«0; k<270; k++){ 

outnode [ il . blc>cK2 [ j ] . layotwgt [k]- 
(doublfc) (randomO % 101)/100.00 - 0.5 - 

) 

) 

) 

for(i-0; 1<4; i++){ /* Normalize initial values by fan in */ 

for(J-«0; J<4, .1++){ 
for(k-0; k<270; k++){ 
outnode [ i] .block2 (j ] layotwgt [k]»« 

(outnode[i].block2[j].layotwgt[k])/1080.00; 

) 

) 

} 

1ay2wgt-l.0/25.0; /* Constant weight for all 2nd hidden layer nodes */ 

/* Initialize four Gabor function weight arrays for first layer 
(maximum value of each function is 0,5 to aid sigmoid computation) */ 

fx-0.0; 
fy-1.0; 
gabcount-0; 
for(1-0; i<-7; i++){ 
for(j-0; j<-7; j++){ 
a- 0.2*(float)i - 0.7; 
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b- 0.2*(£loat)j - 0.7; 
gabor[0)t gabcount]- 

0.5*exp(-pi*(a*a/3.0+b*b/3.0))*cos(2.0*pi*(fx*a+fy*b)) 
gabco'ant»gabcoimt+l; 

) 


fx-1.0 ; 
fy-0.0; 
gabcount-O; 
for(1-0; i<-7; i++){ 
for(j-0; j<-7: j++){ 
a- 0.2*(float)i - 0.7; 
b- 0.2*(float)j - 0.7; 
gabor[1J[gabcount]« 

0.5*exp(-pi*(a*a/3.0+b*b/3.0))*cos(2.0*pl*(fx*a+fy*b)) 
gabcount»"gabcount+l; 

) 


fx=l.4142136; 
fy-1.4142136; 
gabcount-0; 
for(i-0; i<-7; i++){ 
for(j-0; j<-7; j++){ 
a- 0.2*(float)i - 0.7; 
b- C.2*(float)j - 0.7; 
gabor [ 2 J [ gabcount J «• 

0.5*exp( -pl*(.a*a/3.0+b*b/3.0) )*cos (2.0*pi*(fx*a+fy*b)) 
gabcoi*nt“gabcoiint+l; 

) 


fx-1.4142136; 
fy-0.0-1.4142138; 
gabcount-O; 
for(i-0; i<-7; !++){ 
for(j-0; j<-7; j++){ 
a- 0.2*(float)i - 0.7; 
b- 0.2*(float)j - 0.7; 
gabor [ 3 ] [ gabcount ] »• 

0.5*exp(-pi*(a'*a/3.0+b*b/3.0))*cos(2.0*pi*(fx*a4fy’'"b)) 
gabcount»gabcount+l; 

) 


/* End of initialization; begin loops for training and testing. */ 
for(bigloop-1; bigloop<-numruns; bigloop++){ 

/* Begin training loop */ 
for(tmloop-0; tmloop-Cnumtm; tmloop++){ 
randnum-CiandomO % numexvec); 
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/* Begin propagation thru network */ 

/* Begin receptive field window shifting loop -•/ 
laylposn-0; 

for (rows tplx“0; rovstpix<»6656; rowstpix“rowstpix+5.l2) f 

for(colstpix-rowstpixi colstpix<«rowstpix+120; colstpix-colEtpix+4){ 
index-0; 

for(i«»colstpix: l<-=olstpix+896; l-i+128)( 
for(j-i; j<-(i+7), j++){ 

rfarraylindex]-(double)x[(randr.\im*806A)+J] ; 
index-index+1; 

) 

) 

/* LAMBERTIZATION AND CONTRAST NORMALIZATION MODULE (processes each 
receptive field window) */ 
abstun-0.0; 

for(i-0; i<84: i++){ 
absum-absum+rfarray[i]; 

> 

avbrlght-abstuii/6A. 0; 
for(i-0: i<64; i++){ 
rfarraytij“rfarray[i]-avbright; 

) 

es\im-0.0; 

for(i-0; i<64: i++){ 
esuin-esuin+(rfarray [ i]*rfarray[ i]); 

) 

energy»sqrt(esum); 

if (energy !- 0.0) { /* Do not normalize if ei*ergy-0,0 

for(i-0; i<64; i++){ 
rfarray[i]-rfarray[1]/energy; 

) 

) /* end noxnnalization module */ 

/* Propagate thru to first hidden layer outputs */ 

for(i“0; 1<4; i++) { /* Load each 1st layer block from single window 

stim-0.0; 

for(j-0; j<64; j++){ 

sum - sum + gabor[l](j)*rfarray[j]; 

) /* Finished calculating one block, one position (first layer) */ 

if (s’jm<-20.0) {laylout [ i] [laylposn]“0.0;) 

else if(sum>20,0) (laylout[i][laylposnJ-1 0;) 

else (laylout[1][laylposn]»-(double)(1.0/(1.0+exp(-sum)));) 

) /* Finished calculating four blocks, one position (first Icyer) */ 

laylposn-laylposn+l; 

) 

) /* end of receptive field window shifting loops (all positions) */ 
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/* Propagate thru second layer */ 
for(bshlft-0; bshift<4; bshlft++){ 
lay2posn-0; 

for'rowsthl—O; rowsthl<“279; rowsthl»rowsthl+31)( 

for(colsthl“rowsthl; colsthl<“rowsth?.+26; colsth r+){ 
index-0; 

for(i-colsthl; i<-cc.'sthl-«'124; 1-1+31) { 
for(j-l; J<-(i+4); j++){ 
pharray[index]"laylout[bshift][j]; 
i <idex-lndex+l; 

) 

) 

sum-O.0; 

for(1-0; i<25: !++){ 
s\iin-S'uiii+pharray [ i ] ; 

) 

lay2out[bEhtft] [lay2posn] - svim*lay2wgt; 
lay2posn - lay2posn + 1; 

) 

} 

} 

/* Propagate thru output layer */ 
for(1-0; 1<4; !++)( 
sum-O.0; 

for(j-0; j<4: j++){ 
for(k-0; k<270; k++){ 

sum - sum + outnode[i].block2(j).layotwgt[k)*lay2out[j][k]; 

} 

) 

if ((sum+outtheta[i])<-20.0) (y[i]-0.0;) 

else if((sum+outtheta(i])>20.0) {y[i]-1.0;) 

else {y[i]-(double)<1.0/(1.0+exp(-<sum+outtheta[i]))));) 

) /* End calculation for all FOUR output nodes y(4] */ 

/* End propagation thru network */ 

/* Backpropagation to update output weights */ 

for(n-0; n<4; n++){ 

if (class[randnuffl] — n) (d[n]-1.0;) 
else {d[n]-0.0;) 

) /* Filled desired output values, d[4], based on given class */ 

for(n-0; n<4; n++){ 

errorout[n] - y[n]*(l.0-y[n])*(d[n]-y[n]); 

) 

for(i-0: i<4; 1++){ 
for(1-0; J<4: j++){ 
for(k-0; k<270; k++){ 
outnode[l].block2[j].layotwgt(k] - 

outnode(1].blockZfj].li otwgt[k] + 

(eta/1080.0)*errorout[i]*lay2out[j][k] 
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) 

} 

) 

for(n“0; n<^; n++)( 

outtheta[n] - outtheta[n] + (eta/1080.O)*errorout[n]; 

) 

/* end of backprop */ 

) /* end of training loop (tmloop) */ 

/* Begin test loop to run thru each test vector */ 
numrlght-O; 

for(testvec—nrunexvec; testvec<(numexvec+nuinttvec) ; test’'^ec++) { 

/( uegin propagation thru network */ 

/* Begin receptive field window shifting loop */ 
laylposn-0; 

for(rowstplx-0; rowstplx<-6656; rowstplx-rowstplx+512){ 

for(colstplx-:'owstplx; colstplx<-rowstpix+120; colstpix-colstpix+A) { 
index-0; 

for(l-colstpix; i<-colstpix+896; 1—i-t-128){ 
for(j-l; j<-(i+7); j++){ 

rfarray [ lndex] = <'double)x[ (testvec*8064)+j ] ; 
lndex*'index+l; 

) 

) 

/* LAMBERTIZATION AND CONTRAST NORMALIZATION MODULE (processes each 
receptive field window) */ 
absuffl-0.0; 

for(1-0; 1<'A4; 1++) ( 
absum-absum+rfarrayl1]; 

) 

avbright-absuin/64.0; 
for(l‘-0; 1<64; i++) { 
rfarray[1]-rfarray[1)-avbright; 

) 

esum-O.0; 

for(i“0; 1<64; 1++)( 
esuiii»esuic+(rfarray [ i]*rfarray [ 1] ) ; 

) 

energy-sqri(esum); 

if (energy !- 0.0) { /* Do not normalize if energy-0 0 */ 

for(l-0; 1<64; i++){ 
r€array[l]-rfarray[i]/energy; 

) 

) /* end nomallratlon module */ 
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/* Propagate thru to first hidden layer outputs */ 

for(i-0; i<4; i++){ /* Load each 1st layer block from single window */ 

sum-0.0; 

for(J-0; j<64: j++){ 

sum - sum + gabor[i][j]*rfarray[j}; 

) /* Finished calculating one block, one position (first layer) */ 

if (svim<-20.0) {laylout[i] [laylposn]-0.0: } 

else if (su]ji>20.0) (laylout[l] [laylposn]-1.0;) 

else {laylout[i][laylposn]-(double)(1.0/(1.0+exp(-sum)));) 

) /* Finished calculating four blocks, one position (first layer) */ 

laylposn-laylposn+1; 

) 

) /* end of receptive field window shifting loops (all positions) */ 

/* Propagate thru second layer */ 
for(bshlft-0; bshlft<4; bshift++){ 
lay2posn-0; 

for(rowsthl-0; rowsthl<-279; rowsthl-rowsthl+31)( 

for(colsthl-rowsthl; colsthl<-rowsthl+26; colsv:hl++){ 
index-0; 

for(i-colsthl; i<-colsthl+124; i-i+31){ 
for(j-l: j<-(i+4): j++)( 
pharray[lndex]-laylout[bshift][j]I 
index-lndex+1; 

) 

) 

sum-0.0; 

for(i-0; i<25: !++){ 
sum-sum+pharray[1]; 

} 

lay2out[bshift][lay2posn] - sum*lay2wgt; 
lay2posn - lay2posn + 1; 

) 

) 

) 

/* Propagate thru output layer */ 
for(1-0; i<4; !++)( 
sum-0.0; 

for(j-0; j<4; j++){ 
for(k«0; k<270; k++)( 

sum - sum + outnode[1].block2[j}.layotwgt[k]*lay2out[j][k]; 

) 

) 

if((sum+outtheta[l])<-20.r/) {yIl]-0 0;) 

else if((sum+outtheta[i])>20.0) {y[i]-1.0;) 

else {y[i]-(double)(1.0/(1.0+exp(-(sum+outtheta[l]))));) 

} /* End calculation for all FOUR output nodes y[i] */ 

/* End propagation thru network */ 
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/* Decision criteria for correct classification */ 
neurtcnt«0; 
for(n-0; n<4; n++){ 

If((y[n]>0.8) && (class[testvec)—n)) {neurtcnt-neurtcnt+1;) 
if((y[n]<0.2) && (class[teEtvec]I^n)) (neurtcnt'-neartcnt+l;) 

) /* End loop for checking all four output neurons */ 

If(n®urtcnt-“4) (numrlght-n^umright+l;) 

) /* end of testing loop (testvec) */ 

itnum bigloop*nunitm; 

accuracy «• (float)ntiinrlght/(float)n\imttvec; 

printf("Training Epoch #: %d Accuracy-%f\n", itnum, accuracy) 

) /* end of bigloop */ 

) /* end of main */ 
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/* Image clascifIcati.'im netj'crk with hlgh-order (second-order) 
network as ovtput layer. Passes DOUBLE precision values to 
NORMALIZING routine for second-order processing. This 
version exceeded machine memory capacity and was scrapped 
before completion. See later program versions. 

Name: IMAGNET2.C */ 

^include <stdio.h> 

#include •diath.h> 

#include <struct.h> 

^define numexvec 44 
#define numttvec 44 

#define numvec 88 /* Specify total number of vectors */ 

#define nximfaat 8064 /* Specify nximber of features per vector */ 

#deflne pi 3.14159265358979 

struct { 

double x.2[4][270]; 
double xquad[4][36583]; 

) vector[88]; 

char *aialloc(); 

malnO ( 

int *flag, *x, *class, gabcount; 

Int i, j, n, blgloop, Lshift, row, qcotmt, block; 

int rowstpix, colstpix, rovrsthl, colsthl, index, laylposn, layZposn 
double rfarray[64], absxim, avbright, esum, energy, lay2wgt; 
double sum, fx, fy, a, b, sumsqdif, s, mean; 
double laylout[4][434], lay2out[4][270], pharray[25]; 
double gabor[4][64]; 


/* Open and read in data file to arrays flag[], class[], x[]l] */ 
FILE *fp: 

fp-fopen("biglmage","r">; 
if(fp—NULL) exit(O); 

flag - (int *) malloc(n\imvec * sixeot(int)) ; 
class “ (int *) malloc(numvec * sizeef(int)); 

X (int *) malloc(numvec * numfeat * slzeof(int)); 

for(1-0; Knumvec; !++){ 
fscanf(fp,"%d",&flag[l]); 
for(J-0; j^umfeat; j++){ 
fscanf (fp, "%d” ,6tx[ 5.*numfeat+J ]) ; 

) 

fscanf (fp, ''%d" ,6iclass[ i] ) ; 

) 

fclose(fp); 
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/* Test printout module: */ 

/* for (1-0; Knumvec; !++){ 

prIntf("Flag-%d Class-%d\n",flag[1],classI1]); 

) */ 

/* Establish constant "weight” values for Gabor function orientations 
and "phase synchronising" stunmatlon layer */ 

lay2wgt-1.0/25.0; /* Constant weight for all 2nd hidden layer nodes */ 

/* Initialize four Gabor function weight arrays for first layer 
(maxinmm value of each function Is 0.5 to aid sigmoid computation) */ 

fx-0.0; 
fy-1.0; 
gabcount-0; 
for(1-0; l<-7; !++){ 
forCj-O; j<..7; J++){ 
a- 0.2*(float)l - 0.7; 
b- 0.2*(float)j • 0.7; 
gabor[0][gabcount]- 

0.5*exp(-pl*(a*a/3.0+b*b/3.0) )*cos(2.0*pi*(fx’^a+fy*b)) ; 
gabcount-gabcount+1; 

1 

) 

f x-1.0 
fy-0.0; 
gabcount-0; 
for(1-0; i<-7; !++){ 
for(j-0; J<-7; j++){ 
a- 0.2*(float)i - 0.7; 
b» 0.2*(float)J - 0.7; 
gabor[1]t gabcount]- 

0.5*exp(-pi*(a*a/3.0+b*b/3.0> )*cos (2.0'*pl*(fx*a4fy*b)) ; 
gabcount-gabcount+l; 

) 

) 

fx-1 4U2136; 
fy-1.4142136; 
gabcount-0; 
for(1-0; l<-7; i++){ 
fcr(j-0; J<-7; j++){ 
a- C.2*(float)i - 0.7; 
b- 0.2*(float)j - 0.7; 
gabor 12][gabcount]- 

0.5*exp(-pl*(a*a/3.0+b*b/3.0))*cos(2.0*pl*(fx*a+fy*b)); 
gabcount-gabco\mt+l; 

) 

) 
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fx-1.4142136; 
fy-0.0-1.4142136; 
gftbcount-O; 
for(1-0; 1<«7; !++){ 
for(j-0; J<-7; J++){ 

«- 0.2*(float)i - 0.7- 
b- 0.2*(float)j - 0.7; 
gfibor [ 3 ] [ gabcotint ] - 

0. S’^expC -pl*(a*a/3.0+b*b/3.0) )*cos (2.0*pt*(f x*a+fy*b) ) ; 
gabcotint'^gabcount+1; 

} 

) 

/* End of module establisihlng weights. Begin loop thini all vectors. */ 

for(blgloop'-0; bigloopCnumvec; bigloop++){ 

/* Begin propagation thru network */ 

/* Begin receptive field window shifting loop */ 
laylposn«*0; 

for(rowstpix-0; rowstplx<-6656 ; rowstpix-rowstplx-(-512> ( 

for(colstpix-rowstpix; colstplx<-rowstpix+120; colstpix-colstpix+4){ 
index-0; 

for(l-colstplx; l<-colstplx+896; i-i+128)( 
for(j-i; J<-(i+7); J++)( 

rfarray[index)-(doub.i.e)x[ (blgloop*8064)+j); 
index-index+l; 

) 

) 

/* LAMBERTIZATION AND CONTRAST NORMALIZATION MODULE (processes each 
receptive field window) */ 
absum-O.0; 

for(i-0; 1<64; i++){ 
absvun-absum-t-rfarray [ 1] ; 

} 

avbright-ahsum/64.0; 
for(i-0; 1<64; i++){ 
rfarray[i]-rfarray(i]-avbrlght; 

) 

esum-0.0; 

for(1-0; 1<64; !++){ 
esuin-esum-f(rfarray[i]*rfarray[i]); 

) 

energy-sqrt(e8um>; 

if (energy f- 0.0) ( /* Do not normalise if energ/-0.0 */ 

for(1-0; 1<64; i-H-){ 
rfarray[1]-rfarray[i]/energy; 

) 

) /* end normalization module */ 
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/* Propagate thru to first hidden layer outputs */ 


for(i"iO; 1<4: i++){ /* Load each 1st layer block from single window */ 

sum-0.0; 

for<j-0; J<64: 

sum - sum + gabor[l][j]*rfarray[JJ; 

) /* Finished calculating one block, one position (first layer) */ 

if(sum<-20.0) {laylout[i][ldylposn]-0.0;) 

else if (s\im>20.0) (layloutli] [laylposn]-1.0:) 

else (laylout[i][laylposn]-(double)(l.0/(1.0+exp(-sum)));) 

) /* Finished calculating four blocks, one position (first layer) */ 

laylposn-laylposn-t-1; 

) 

) /* end of receptive field window shifting loops (all positions) */ 

/* Propagate thru second layer */ 
for(bshift-0; bshift<4; bshift++){ 
lay2posn-0; 

for(rowsthl"0; rowsthl<-279; rowsthl-rowsthl+31){ 

for(colsthl-rowsthl; colsthl<-rowsthl+26; colsthl++){ 
index-0; 

for(l-colsthl; l<-colsthl+124; l-i+31){ 
for(j-i; j<-<i+4); j++){ 
pharray[index]-laylout[bshift][J]; 
index-index-f 1; 

) 

) 

sum-0.0; 

for(1-0; 1<25; !++){ 

Evun-sum+pharray[ 1] ; 

) 

lay2out[bshift][lay2posn] - sum*lay2wgt; 
lsy2posn - lay2posn + 1; 

) 

) 

) 

/* End propagation of one vector thru output of summation layer */ 

/* For each run thru blgloop (over each new vector), fill structure 
array vectorI88].x2[4][270] with propagated results from outputs 
of summation layer. The filled x2 array will be the input for 
the final, high-order layer. */ 

for(i-0; i<4; !++){ 
for(J-0; J<270; J++)( 

vector[blgloop).x2[l]tj] - lay2out[i][j]; 

) 

) 

) /* end of blgloop */ 
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/* Begin hlgh-order modules */ 

/* Fill arrays for second-order Input combos 
vector[88].xquad[A3[36585] */ 

for(rowO; row<numvec; row++){ 
for(block“0; block<4; block++){ 
qcount->0; 

for(1-0; 1<270; !++)( 
for(J-l; j<270; j++){ 
vector[rov].zquad[block][qcount] — 

vector[row].x2[block][1] * vector[row].x2[block][j] 
qcount-qcount+l; 

) 

) 

) 

) 


) /* end of main */ 
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/* Image classification network with high-order (second-order) 
network as output layer. Passes DOUBLE precision values 
for second-order processing. Second-order combinations 
are mul'tiplied among the 270 nodes within each of the 
four separate orientation groups of the phase 
synchronization layer outputs. 

Name: IMAGNET3.C */ 

#include <stdlo.h> 

#lnclude <math.h> 

#include <struct.h> 

#define numexvec 44 
#define niunttvec 44 

#deflne numvec 88 /* Specify total number of vectors */ 

#deflne numfeat 8064 /* Specify ntimber of features per vector */ 

#deflne n\imruns 60 /* Specify number of result plot points */ 

#define numtm 1000 /* Specify # training iters between test/plot */ 

#define eta 0.35 

#define pi 3.14159265358979 

struct { 

double x2[4][270]; 

) vectorI88]; 

struct ( 

double linwgt[4][270]; 
double quadvgt[4][36585] ; 

) neuron[4]; 

char *malloc<); 

mainO { 

int *flag, *x, *class, gabcount. fulloop, tmloop, randnum, neurtcnt; 
Int i, j, n, bigloop, bshift, row, qcount, block, testvec, numright; 
int rowstpix, colstpix, rovsthl, colsthl, index, laylposn, lay2posn; 
int Itnum; 

double rfarray[64], abs\un, avbright, esum, energy, lay2vgt; 
double sum, fx, fy, a, b, 1 insum, qtiadsum, y, d, error; 
double laylout[4][434], lay2out[4][270], pharray[25]; 
double gabor[4][64], theta[4], xquad[4][36585]; 
float accuracy: 


/* Open and read in data file to arrays flag[], class[], x[][] */ 
FILE *fp; 

fp-fopen("biglmage","r"); 
if(fp—NULL) exlt(O); 
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flag "■ (Int *) mallocCnumvet * wireof(Int)); 
class - (int *) inalloc(nujnvec * sizeof (ipt)); 

X -• (Int *) malloc(nvunvec * numfeat * sizeof (Int)); 

for(i»«0; Knumvec; i++){ 
fEcanf(fp,*‘%d",&flag[i]); 
for(j-0; J<numfeat; J++>{ 
fscanf (fp, "%d" ,&xti'^vunfeat+j ] ) ; 

) 

fscanf(fp,"%d",&cl8ss[l]); 

) 

fclose(fp); 

/* Test printout module; */ 

/* for(i-0; Knumvcc; !++){ 

printf("Flag“%d Class«%d\n",flag[i],class[i]); 

) */ 

/* Establish constant "weight" values for Gabor function orientations 
and "phase synchronizing" svimmation layer */ 

lay2wgt«-l.C/25.0; /* Constant weight for all 2nd hidden layer nodes */ 

/* Initialize four Gabor function weight arrays for first layer 
(mAxlmum value of each function is 0.5 to aid sigmoid computation) */ 

fx-0.0; 
fy-1.0; 
gabcount«0; 
for(i-0; l<-7; i++)( 
for(j-0; j++){ 

a- 0.2*(float)i - 0.7; 
b- 0.2*(float)j - 0.7; 
gabor[0][gabcount]« 

0.5*exp(-pi*(a*a/3.0+b*b/3.0))*cos(2.0*pi*(fx*a+fy*b)); 
gabcount-gabcount+1; 

) 

) 

fx-l.0; 
fy-0.0; 
gabcount“0; 
for(1-0; l<-7; !++){ 
for(J-0; j<-7. j++){ 
a- 0.2*(float)l - 0.7; 
b- 0.2*(float)J - 0.7; 
gabor[1](gabcount]- 

0.5*exp(-pl*(a*a/3.0+b*b/3.0))*cos(2.0*pi*(fx*a+fy*t)), 
gabcount-gabcount+1; 

) 

) 
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fx-1.4U2136; 
fy-l.4142136; 
gabcounti-0; 
for(1-0; l<-7; !++){ 
for(j-0; j<-7; J+t){ 
a- 0.2*(float)i - 0.7; 
b- 0.2*(float)J - 0.7; 
gabor[2][gabcount]- 

0.5*exp('pl*(a*a/3.0+b*b/3.0))*cos(2.0*pi*(fx*a+fy*b)); 
gabcoTint-^abcowt+l; 

) 

) 

fx-l.4142136; 
fy-0.0-1.4142136; 
gabcount-0; 
for(i-0; l<-7; !++){ 
for(j-0; j<-7; j++){ 
a- 0.2*(float)i - 0.7; 
b- 0.2*(float)j - 0.7; 
gabor[3][gabcount]- 

0.5*sxp(-pi>(a*a/3.0+b*b/3.0))*cos(2.0*pi*(fx*a+fy*b)); 
gffibcount-gabcotuit+1; 

) 

) 

/* End of module establishing weights. Begin loop thru all vector.*;. */ 

for(bigloop-0; bigioop<numvec; b5gloop++){ 

/* begin propagation thru network */ 

/* Begin receptive field window shifting loop */ 
laylposn-0; 

for(rowstpix»0; rowstplx<»6656; rowstptx-rowstpix+512){ 

for<colstpix-rowstpix; colstpix<-rowstpix+120; colstplx-colstpix+4){ 
index-0; 

for(i-colstpix; i<«coIstpix+896; i-i+128){ 
for(J-l; J<"(i+7); j++M 

rfarray[index]-(double)x[<blgloop*8064)+j); 
index-index-t-l; 

) 

) 

/★ LAMBERTIZATION AND CONTRAST NORMALIZATION MODULE (processes each 
receptive field window) */ 
abs\UB-0.0; 

for(1-0; 1<64; i++){ 
absun-abstUQ4rfarray[i]; 

) 

avbright-absum/64.0; 
for(i-0; 1<64; !++){ 
rfarray[l]-rfarray[i]-avbright; 

) 
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esum^O.0; 

for(1-0; 1<64; t++){ 
esuin»*esuni+(rfartay[l]*rfarrav[l]); 

) 

energy-isqrt(esuni); 

If (energy I- 0.0) ( ** Do not Aormalize If energy-0.0 */ 

for(1-0; l<a4; H+) { 
rfarray f1j-rfarray[1]/energy; 

) 

) /* end normalization module */ 

/* Propagate thru to first hidden layer outputs */ 

for(i-0; 1<4; 1++){ /* Load each 1st jiayer block from single window */ 

svuc-a. 0; 

for(j»0; j<64; j++)( 

sum • sum + gaborf1][j]*rfarray[j]; 

) /* Finished calculating one block, one position (first layer) */ 

if(5um<-20.0) {layloucti][laylposn]-0.0;) 

else if (s\un>?.0.0) {laylout [ 1] [laylposnj-l. 0;) 

e.’^e {laylouttl] [laylposn]‘ (double)(1.0/(1.0Hei£p(-stm))):) 

) /* Finished calculating four blocks, one position (first layer) */ 

laylposn-laylposn+l; 

J 

) /* end of receptive field window shifting loops (all positions) */ 

/* Propagate thru second layer */ 
for(bshift-0; bshift<4; bshif_++){ 
lsy2posn-0; 

for(rowsthl“0; rowsthl<-279; rowsthl-rowsthl+31){ 

for(colsthl-rowsthl; colsthl<-rowsthl+26; colsthl++>{ 
index-0; 

for(l-colsthl; i<-colsthl+124; i-l+31)l 
for(j-l; j<-(i+4): j++){ 
pharray[index]-laylout[bshlft](j]; 
lndex-index-t-1; 

) 

) 

sum-0.0; 

for(l-0: 1<25; !++){ 
sum-sum+pharray[1]; 

) 

lay2out[bshift][lay2posn] » sum*lay2wgt; 
lay2posn - lay2posn + 1; 

) 

) 

) 

/* End propagation of one vector thru output of summation layer */ 
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/* For each run thru blglpop (over each new vector), fill structure 
array vector[88].x2[A][270] with propagated results from outputs 
of stumatlon layer. The filled x2 array will be the input for 
the final, high-order layer. */ 

for(i»0; 1<4: i++){ 
for(j-0; j<270; j++){ 

vector[bigloop].x2[i][j] •• lay2out[l][j]; 

) 

) /* end of bigloop */ 

/* Begin hlgh-ordar modules */ 

/* Initialization */ 

srandomC(unsigned)time(NULL)); /* Inlt randomO seed off clock */ 

for(n-0; n<4; n++){ 

forCblock-0; block<4: block++){ 
for(j-0; j<270; j++){ 

neuron[n].llnwgt[block][j) » (double)(randomO % 101)/100.00 - 0.5 
neuron[n].linwgt[block][j] - neuronfn].linwgt[block][i]/147420.00; 
} 

) 

. ) 

for(n-0; n<4; n++){ 

for(block-0; block<4; block++){ 
for(j-0; j<36585; j++){ 

neuron[n].quadwgt[block][j]-(double)(randomO % 101)/100.C0 - 0.5; 
neuron[n].quadwgt[block][j]-neuron[n].quadwgt[block][j]/147420.00; 
) 

) 

} 

for(n-0; n<4; n++)( 

theta[n] - (double)(randomO % 101)/100.00 - 0.5; 
theta[n] - theta[n]/147420.00; 

) 

/* End of initialization; Begin training and testing loops */ 

for(fulloop-1; fulloop<-numruns; fulloop++){ 
for(tmloop-0; tmloop<numtm; tmloop-M-)( 

randnum-(randomO t nxunexvec); /* Select random exemplar vector */ 

/* Multiply out second-order combinations for present exemplar */ 
for(block-0; block<4; block++){ 
qcount-0; 

for(i-0; 1<270; !++){ 
for(J-i; J<270; J++)( 
xquad[block][qcount] » 

vjjctor[randn\im] .x2 [block] [1] * vector[randnum] .x2 [block] ( j ) 
qcovint-qcount+l; 

) 

) 

) 
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for(n-0; n<4; n i-+) { /* Loop to change neuron */ 

1lnsum-0.0; 

forCblock»0: block<4; block++){ 
tor(j>-0; j<270; j++){ 
lins'jm - llnsm + 

neuron[n].llnwgt[block][j]*vfictor[randntm].x2[block][j]; 

) 

) /* end block */ 
quadsvuii»0.0; 

for(block-0; block<4: block++){ 
for(j-0; j<36585; ]++){ 

quadsum - quadsum + neuron[nl.quadwgt[block][j]*xquad[block][j]; 

) 

) /* end of block */ 

sun > llnsum -i- quadsum; 

if ((sum+theta[n])<-20.0) (y-O.O;) 

else If ((sum+theta[n])>20.0) (y-l.O;) 

else (y-(double)(1.0/(1.0+exp(-(sum+theta[n]))));) 

/* Criteria for training each class to fire a specific neuron +/ 
if (class[randniim]—n) (d-i.O;) 
else (d-O.O;} 

error-y*(1.0-y)*(d-y); 

/* Update weights */ 
for(block-0; block<4; block++){ 
for(J-0; j<270; j++){ 

neuron[n].llnwgt[block][j] - neuron[nl.llnwgt[block](j] + 

(eta/147420.0C)*error*vector [randniim] . x2 [block] [ j ] ; 

) 

) 

for(block-0; block<4; block++){ 
for(j-0; j<36585; J++){ 

neuron[n].quadwgt[block][J] - neuron[n].quadwgt[block][j] + 

(eta/la7420.00)*error*xquad[block][j]; 

) 

) 

theta[n] - theta[n] + (eta/147420.00)*error; 

) /* End of neuron (n) changing loop */ 

] /* End of tmloop */ 

/* Begin test loop to run thru each test vector */ 
munrlghr.-O; 

for(testvec—numexvec; testvec-cnumvec; testvcc++){ 
neurtcnt-0; 

/* hultlply our second-order combinations for present test vector */ 
for(block-0: block<4; block++)( 
qcount-0; 

for(l-0; 1<270; !++)( 
for(j-l; J<270: j++){ 
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xquad[block] [qcount] •* 

vector[testvec].x2[block][i] * vfcctor[testvec].x2[block][j] 
qcount“qcount+l; 

) 

) 

) 

for(n-0; n<4; n-f+){ /* Loop to change neuron */ 

llns\iffl~0.0: 

for(block-0; block<4; block++){ 
for(j-0: j<270; j++){ 
linsum - 1insum + 

neuron[n].linwgt[block][J]*vector[testvec].x2[block][j]; 

} 

) /* end block */ 
quadsum-0.0; 

for(block“0; block<4; block++){ 
for(j-0; j<36585; j++){ 

quadstiffi - quadsum + neuron[n].quadwgt[block][j]*xquad[block][j]; 

) 

) /* end of block */ 

suw » linsum + quadsum; 

if ((sum+theta[n])<-20.0) {y«0.0;) 

else if ((sum+theta[n])>20.0) (y-l.O;) 

else {y«(double)(1.0/(1.0+exp(-(sum+thcta[n]))));) 

/* Test decision criteria for identifying classes */ 
if ((y>0.8) && (class[testvec]»-n)) {neurtcnt-neurtcnt+1;) 
if ((y<0.2) && (class[testvec]!-n)) {neurtcnt»neurtcnt+l; ) 

) /* End of netiron (n) changing loop */ 
if (neurtcnt“4) (numright-numrlght+l;) 

) /* End of test loop running thru each test vector (testvec) */ 
itnum-fulloop*numtrn; 

accuracy-(float)numright/(float)numttvec; 

prlntf("Training Epoch #; %d Accuracy-%f\n", itnum, accuracy); 

) /* End of fulloop */ 

) /* end of main */ 
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/* Image classification niitwork with high-order (second-order) 
network as output layer. Passes DOUBLE precision values 
for second-order processing. Multiplies out second-order 
combos among the four orientations for each of the 270 
nodes (per orientation) feeding final layer. 

Name: IMAGNET4.C */ 

#lnclude <stdio.h> 

#lnclude '^ath.h> 

#lnclude <stj:uct.h> 

♦define numexvec 44 
♦define numttvec 44 

♦define numvec 88 /* Specify total number of vectors */ 

♦define nvimfeat 8064 /* Specify number of features per vector */ 

♦define numruns 1000 /* Specify number of result plot points */ 

♦define nTirntm 1000 /* Specify ♦ training Iters between test/plot */ 

♦define eta 1,0 

♦define pi 3.14159265358979 

struct { 

double x2[4][270]; 

) vector[88]; 

struct { 

double linwgt[4][270]; 
double quadwgt[270][10]; 

} neuron[4]; 

char *malloc(); 

malnO ( 

Int *flag, *x, *class, gabcount, fulloop, tmloop, randnum, neurtcnt; 
int 1, J, n, blgloop, bshlft, row, qcount, block, testvec, numright; 
int rowstpix, colstplx, rowsthl, colsthl, index, laylposn, lay2posn: 
int itnum, blockl, blockj, spot; 

double rfarray[64], absum, avbright, esum, energy, lay2wgt; 
double sum, fx, fy, a, b, llnstim, q\iadsum, y, d, error; 
double laylout[4][434], lay2out[4][270], pharray[25]; 
double gabor[4][64], theta[4], xquad[270][10]; 
float accuracy; 


/* Open and read In data file to arrays flag[], cla5s[], x[][] */ 
FILE *fp; 

fp-fopen("bigimagf;", "r") ; 
if(fp«-NULL) exlt(O): 
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flag “ (Int *) mallocCnumvec * slzeof(Int)); 
class « (Int *) malloc(ntunvec * slzeof(int)); 

X (int *) malloc(numvec * numfeat * slzeof (int)); 

for(l«0; Knviaivee; 1++) ( 
fscanf(fp,"%d",&flag[i]); 
for(j-0; j<nuiafeat; j++){ 
fscanf(fp,"%d",&x[i*n\uofeat+j]); 

) 

fscanf(ft.."%d".&class[l]); 

) 

fclose(fp); 

/* Test printout module: */ 

/* for(i“0; Knumvec; !++){ 

printf('’Flag»%d Class«%d\n",flag[i] ,class[i]) ; 

) */ 

/* Establish constant "weight" values for Gabor function orientations 
and "phase synchronizing" summation layer */ 

lay2wgt“1.0/25.0; /* Constant weight for all 2nd hidden layer nodes 

/* Initialize four Gabor function weight arrays for first layer 
(maximum value of each function is 0.5 to aid sigmoid computation) */ 


fx-0.0; 
fy-1.0; 
gabcount-0; 
for(i-0; i<-7; i++){ 
for(j-0; j<-7; j++)( 
a- 0.2*(float)i - 0.7; 
b- 0.2*(float)j - 0.7; 
gabor[0][gabcount]- 

0.5*exp(-pi*(e*a/3.0+b*b/3.0))*cos(2.0*pi*(fx*a+fy*b)); 
gabco\snt-gabcount+l; 

) 

) 


fx-1.0: 
fy-O.0 ; 
gabcount“0; 
for(i-0; i<-7; i-^« ){ 
for(J-0; J<-7; ,|++){ 
a- 0.2'^(float)l - C.?; 
b- 0.2*(float)j - 0.7; 
gabor[1][gabcount]- 

0.5*e:xp(-pi*(a*a/3.0+b*b/3.0))*cos(2.0*pl*(fx*a+fy*b) ) ; 
gabcount-gabcount+l; 

) 


2 






l'x-1.4142136; 

£y-l.4142136; 
gabcovmtoO; 
for<i-0: i<-7; !++){ 
for(J-0; j<-7; 

•- 0.2*(flo«t)l - 0.7; 
b- 0.2*(float)j - 0.7; 
gabor [ 2 ] [ gabcount ] 

0.5*exp(-pi*<a*a/3.0+b*b/3.0))*cos(2.0*pl*(fx*a+fy*b)); 
gabcount-xgabcount^l; 

} 

) 

fx-1.4142136; 
fy-0.0«1.4142136; 
gabcount-O; 
for(i-0; l<-7; !++){ 
for(j-0; j<-7; j++){ 
a- 0.2*(float)! - 0.7; 
b- 0.2*(float)J - 0.7; 
gabor [ 3 ] [ gabcourit ] • 

0.5*exp<-pi*<a*a/3.0+b*b/3.0))*cos(2.0*pi*(fx*a+fy*b)); 
gabcount“gabcount+l; 

) 

) 

/* End of module establishing weights. Begin loop thru all vectors. */ 

for (bigloop-0; bigloop-Cnumvec; bigloop++) { 

/* Begin propagation thru network */ 

/* Begin receptive field window shifting loop */ 
laylposn-0; 

for(rowstpix-0; rowstpix<-6656; rowstpix™rowstpix;-512) { 

for(colstplx-rowstpix; colstplx<-rowstplx+120; colstpix-colstpix+4){ 
index>-0; 

for<i“COlstpix; i<«-colstpix+896; i-i+128){ 
for(j-l; j<-(i+7); j++){ 

rfarray[index1-(double)x[(blgloop»^8064)+j]; 
lndex»index-i-l; 

) 

) 

/* LAMBERTIZATION AND CONTRAST NORMALIZATION MODULE (processes each 
receptive field window) */ 
absumoO.0; 

for<1-0; i<64; 14+)( 
absum-absun+rfarray(1]; 

) 

avbrlght-absum/64.C; 
for(i-0; 1<64; !++)( 
rfarray[1]-rfa rray[i]-avbrIght; 

} 
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esum-O.0; 

for(1-0; 1<64: !++)( 
estun-«isum-f(rfarray [l]*r£array[l]) ; 

) 

energy-sqrt(esu»); 

If (energy I- 0.0) { /* Do not normalize If energy-0.0 */ 

for(1-0; 1<64; !++)( 
rfarray[1]-rfarray[1]/energy; 

} 

) /* end nomnallzatlon module */ 

/* Propagate thru to first hidden layer outputs */ 

for(i-0; 1<4; !++){ /* Load each 1st layer block from single window */ 

sum-O.0; 

for(j-0; J<64; J++){ 

sum - sum + gabor[l][j]*rfarray[J]; 

) /* Finished calculating one block, one position (first layer) */ 

if(sum<-20,0) (laylout[l][laylposnj-0.0;) 

else lf(sum>20.0) {laylout[ll[laylposn]-1.0;) 

else (laylout[i][laylposn]-(double)(1.0/(1.0+exp(-sum)));) 

) /* Finished calculating four blocks, one position (first layer) */ 

laylposn-laylposn-f 1; 

) 

) /* end of receptive field window shifting loops (all positions) */ 

/* Propagate thru second layer ^/ 
for(bshtft-0; bshift<4; bshi£t++){ 
l£y2posn-0; 

for(rowsthl-0; rowsthl<-279; rowsthl-rowsthl+31){ 

for(colsthl»rowsthl; colsthl<-rowsthl+26; colsthl++){ 
index-0; 

for(i-colsthl; K-co]sthl+124; i-i+31)( 
for(j-i; J<-(i+4); j++)( 
pharray[index]-laylout[bshift][j]; 
lndez-lndex-i-1; 

} 

) 

axun-0.0; 

for(1-0; 1<25; !++){ 

8uiD-sum-»-phari.'ay[i]; 

) 

lay2out[bshift][lay2posn] - 8um*lay2wgt; 
lay2posn - laylposn -f 1; 

) 

) 

) 

/* End propagation of one vector thru output of summation layer */ 
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/* For «ach rtui thru blgloop (over each new vector), fill structure 
array vector[88].x2[4][270] with propagated results from outputs 
of summation layer. The filled x2 arrey will be the Input for 
the final, hlgh-order layer. */ 

for(1-0; 1<4: 1++){ 
for(j-0; J<270; 

vector[blgloopl.x2[l][j] - lay2out[ll[J]; 

) 

) 

) /* end of blgloop */ 

/* Begin hlgh-order modules */ 

/* Initialization */ 

srandom((unsigned)time(NULL)); /* Inlt randomO seed off clock */ 

for(n-0; n<4; n++){ 

for(block-0; block<4; block++)( 
for(j-0; J<270; J++)( 

neuron[n].linwgt[block][j] - (double)(random() % 101)/100.00 - 0.5 
neuron[n].llnwgt[block][j] - neuron[n].llnwgt[block][j]/147420.00; 
) 

) 

) 

for(n«0; n<4; n++){ 

for(spot-0; spot<270; spot++){ 
for(j-0; J<10: j++){ 

neuron[n].quadwgt[spot][j]-(double)(random() % 101)/100.00 - 0.5; 
neuron[n].quadwgt[spot][J]-neuron[n].q\iadwgt[spot][j]/147420.00; 

) 

} 

) 


for(n-0; n<4; n++){ 

theta[n] - (double)(random() % 101)/100.00 - 0.5; 
theta[n] - thftta[n]/147420.00; 

) 

/* End of Initialization; Begin training and testing loops */ 

fox(fulloop-l; fulloop<-numruns; fulloop++)( 
for(tmlcop-0; tmloop<numtm; tmloop++){ 

randnum-(randomO % numexvec) ; /* Select random *^xemplar vector */ 

/* Multiply out second-order combinations lor present exemplar */ 
for(spot-0; 8pot<270; spot++)( 
qcount-0; 

for(blockl-0; blockKA; blockl++>{ 

for(bloc,kj-blockl; blockj<4; blockj4+){ 
xquad [ .spot ] [qcount ] - 

vector[randnum].x2[blockl][spot] * vector[randnam].x2[block]][spot] 
qcount-qcount+1; 

) 

) 

) 
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/* Loop to change neuron */ 


for(n-0; n<4; n++){ 
llnstun^O.O; 

for(block-0; block<4; block++){ 
for(j-0; j<270; J++){ 
llnsTun — I insum + 

neuron[n].llnwgt[block][j]*vector[randnum].x2[block][J]; 

) 

) /* end block */ 
quadsum-0.0; 

for(spot-0; spot<270; spot++)( 
for(j-0; j<10; j++){ 

quadsum - quadsum + neuron[n] .quadvgt[spot] [J]^quad[spot] [ j] ; 

} 

) 

sum - llnsum + quadsum; 

if ((sum+theta[n])<-20.0) (y-O.O;) 

else if ((sum+theta[n])>20.0) (y-l.O;) 

else {y-(double)(1.0/(1.0+exp(*(s\im+theta[n]))));} 

/* Criteria for training each class to fire a specific neuron */ 
if (class[randnum]“-n) (d-l.O;) 
else (d-0.0;) 

error-y*(1.0-y)*(d-y); 

/* Update weights */ 
for (block-0; block<4; block+-+){ 
for(j-0; j<270; j++){ 

neuron[n].llnwgt[block][j] - neuron[n].llnwgt[block][j] + 

(eta/147420.00)*error*vector[randnum].x2[block)[j]; 

} 

) 

for(spot-0; spot<270; spot++){ 
for(j-0; J<10; J++){ 

neuron[n].quadwgt[spot][j] - neuron[n].quadwgt[spot][j] + 

(eta/147420.00)*error*xquadI spot)[j]; 

) 

) 

theta[n] - theta[n] + (eta/147420.00)*error; 

) /* End of neuron (n) changing loop */ 

) /* End of tmloop */ 

/* Begin test loop to run thru each test vector */ 
numrlght-0; 

for(testvec—nomexvec; testvec<numvec; te8tvec++){ 
neurtcnt-0; 

/* Multiply out second-order combinations for present test vector */ 

for(spot-0; spot<270; spot++)( 

qcount-0; 

for(blockl-0; blockl<4 ; blocki+-+){ 

for(blockJ-blockl; blockJ<4; blockJ++)( 
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squad[spot][qco\mt] - 

vector[testvec].s2[blockl][spot] * vector[testvec].x2[biockj][spot] 
qcount-qcount+1; 

) 

) 

) 

for(n“0; n<4; n-H-){ /* Loop to change neuron */ 

linsum-iO. 0; 

£or(blocki-0; block<4; block++)( 
for(J-0; J<270; J++){ 
llnsTJffi llnsum + 

neuron[n].llnwgt[block][j]*vector[testvec].x2[block][j]: 

) 

) /* end block */ 
quadsuin-0.0; 

for(spot-0; spot<270; spot++)( 
for(j-0; j<10; j+-ir){ 

quadsum - quads\ua -t- neuron[n] .quadvgt[spot] [J]*xquad[spot] [ j] ; 

) 

) 

sum - llnsum + quadsum; 

If ((sum+theta[n])<-20.0) (y-O.O;) 

else if ((sum+theta[n])>20.0) (y-l.O;) 

else (y-(double)(1.0/(1.0+exp<-(sum+theta[n]))));) 

/* Test decision criteria for Identifying classes */ 
if ((y>0"8) && (class[testvec]—n)) {neurtcnt-neurtcnt+1;) 
if ((y<0.2) && (class[testvec]!-n)> (neurtcnt-neurtcnt+1;) 

) /* End of neuron (n) changing loop */ 
if (neurtcnt—A) (numright-numright+l;) 

) /* End of test loop running thru each test vector (testvec) */ 
itnum-fulloop*numtm; 

accuracy-(float)numright/(float)numttvec; 

priatf("Training Epoch #: %d Accuracy-%f\n", itnum, accuracy); 

) /* End of fulloop */ 

} /* end of main */ 
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/* Image classification network using multilayer perceptron for final 
output processing; variable number of hidden layer nodes. 

Name: IMAGNET5.C */ 

#include <8tdio.h> 

#include <iiath.h> 

#define numexvec 44 /* Specify number of exemplar vectors */ 

#define numttvec 44 Specify number of test vectors */ 

<!^define n\imfeat 8064 /* Specify nxunber of features per vector */ 

#define midnodes 30 /* Specify number of hidden layer nodes */ 

#define numruns SO /* Specify number of result plot points */ 

#define numtm 1000 /* Specify # training epochs between test/plot */ 

#define eta 1.0 /* Specify learning rate */ 

#define pi 3.14159265358979 

char *malloc(); 

icalnO ( 

Int *flag, *x, *class, neurtcnt, itnum, gabcount, loadcnt; 

int i, J. k, n, bigloop, trnloop, testvec, numrlght, randnum, bshift; 

int rowstpix, colstplx. rowsthl, colsthl, index, laylposn, lay2posn; 

double rfarray[64], absum, avbrlght, esum, energy, lay2wgt: 

double outtheta[4], midtheta[mldnodes], sum, fx, fy, a, b, errsum; 

double laylout[4][434], lay2out[4][270], pharray[25]; 

double gabor[4][64], input[1080], outwgts[4][midnodes]; 

double midwgts[midnodes][1080]: 

double d[4], errorout[4], errormid[midnodes], yout[4], ymid[midnodes] ; 
float accuracy; 

/* Open and read in data file to arrays flag[], class[], x[][] */ 

FILE *fp; 

fp»fopenC"blglmage","r"); 
if(fp—NULL) exlt(O); 

flag - (int *) malloc((numexvec+numttvec) * sireof(Int)); 
class - (int *) malloc((n\ifflexvec-t-numttvec) » slteof (int)) ; 

X » (int *) malloc((numexvec-«-numttvec)*nufflfeat * sizeof (int)) ; 

for (1-0; i<(numexvec-i-numttvec) ; !++)( 
fscanf (fp,"%d" ,«.flag[ 1)) ; 
for(j-0; J<numfeat; J++){ 
fscanf(fp,"%d",&x[i*numfeat+j]); 

) 

fscanf(fp."%d",«iclas8[l]); 

) 

fclose(fp); 

/* Test printout module: */ 

/* for(i“0; l<(numexvec+numttvec); !++)( 
prlntf("Plag-%d\n",flag[l]); 

) V 
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srandomC(unsigned)time(NULL)); /* Init random() seed off clock */ 
/* Initialize output theta and weight arrays */ 
for(i-0: 1<4; i++){ 

outtheta[i] - (double)(random() % 101)/100.00 - 0.5; 
outtheta[l] » outtheta[i]/(double)mldnodes; 

) 

for(i«0; i<nidnodes: 1++){ 

midtheta[i] - (double)(random() % 101)/100.00 - 0.5; 
mldtheta[i] » midtheta[l]/1080.00; 

) 

for(i-0; i<4; !++){ 

for(J-0; J-Cnidnodes; j++)( 

outwgts[l][j] - (double)(random(> % 101)/i00.00 - 0.5; 
outwgts[l][j] - outwgts[1][j]/(double)midnodes; 

) 

) 

for(l«»0; Kmidnodes; 1++) ( 
for(j-0; j<1080; j++){ 

midwgts[i][j] - (double)(randomO % 101)/100.00 - 0.5; 
midwgts[l][j] - midwgtsti]Ij]/1080.00; 

) 

) 


lay2wgt-l.0/25.0; /* Constant weight for all 2nd hiddisn layer nodes ■>’'/ 

/* Initialize four Gabor function weight arrays for first layer 
(maximum value of each function is 0.5 to aid sigmoid computation) */ 

fx-0.0; 
fy-1.0; 
gabcount->0; 
for(i-0; l<-7; !++){ 
for(j-0; j<-7; j4-t){ 
a- 0.2*(float)l - 0.7; 
b- C.2*(float)j - 0.7; 
gabor(0][gabcount]- 

0.5*exp(-pi*(a*a/3.0+b*b/3.0))*cos(2.0*pi*(fx*a+f y*b)); 
gabcountogabcount-tl; 

) 


fX"1.0; 
fy-0.0; 
gabcount»0; 
for(1-0; l<-7; !++)( 
for(j-0; J<-7; j++){ 
a- 0.2*(float)l - 0.7; 
b- 0.2*(float)J - 0.7; 
gabor[1)[gabcount]- 
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0.5*8xp(-pi*(«*a/3.0+b*b/3.0))*cos(2.0*pi*(fx*a+fy*b)); 
gabco\int~gabct'unt-t-l; 

) 

) 

fx-1.4142136; 
fy-1.4142136; 
gabcount-O; 
for(1-0; l<-7; !++){ 
for(j-0; j<-7; j4+){ 
a- 0.2*(float)i - 0.7; 
b- 0.2*(float)j - 0.7; 
gabor[2][gabcount]- 

0.5*exp(-pl*<a*a/3.0+b*b/3.0))*co.-5(2.0*pl*(fx*a+fy*b)) ; 
gabcoimt«gabco\int+l; 

) 

) 

fx-I.4142136; 
fy-0,0-1.4142136; 
gabcount-0; 
for(1-0; l<-7; !++){ 
for(J-0; j<-7; j++){ 
a- 0.2*(float)l - 0.7; 
b- C.2*(fioat)j - 0.7; 
gabor F 3][gabcount]- 

0.5*exp(-pi*(a*a/3.O+b^b/3.0>)*cos(2.0*pi*(fx*a+fy*b)); 
gabcoTint-gabcount-f 1; 

) 


/* End of initialization; begin loops for training and testing. 
for(bigloop-l; blgioop<--numrur.s ; bigloop++)( 

/* Begin training loop */ 

for(tmloop-0; tmloopoiumtim; trnloop++){ 

randnu!n“(random() % numexvec); 

/* Begin propagation thru network */ 

/* Begin receptive field window shifting loop ♦/ 
laylposn-O; 

for(rowstplx“0; rowstpix<-66.S6 ; rowstt)lx-rowst:pix+;'12) { 

for(colstplx-towstpix; col.stplx<«rowstplx+120; colstpi::-colstpij;+4) { 
lndex-0; 

for(i-colstplx; l<-colstpix+95b; i-l+128){ 
for(J-l; j<-(i+7); j++){ 

rfarray [ Index]' (double)x[ (r8ndnuai*8064)+J ] ; 
lndex-lndex+1; 

) 

} 
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/* LAMBERTIZATION AND CONTRAST NORMALIZATION MODULE (processes each 
receptive field window) */ 
absum^O.0; 

for(1-0; 1<64; !++){ 
abs'un-absTim-frfarrayEl]; 

) 

8vbright--absTun/64.0; 
for(1-0; t<64; 1++){ 
rfarray[l]-rfarray[i]-avbrlght; 

) 


estuc-O.O; 

for(l-0; i<64; l-M-){ 
esiijn-esum+(rfarrayti]*rfarray[l]); 

) 

energy-sort(esum ); 

if (energy I- 0.0) { /* Do not normalize if energy-0.0 */ 

for(1-0; 1<64; !++){ 
rfarray[1]-rfarray(1]/energy; 

) 

) /* end normalization module */ 

/* Propagate thru to first hidden layer outputs */ 

for(l-0; 1<4; i++){ /* Load each 1st layer block from single window */ 

sum-O.0; 

for(j-0; 3<64; j++){ 

sum - sum + gaborf1][jl*rfarray[j]; 

) /* Finished calculating one block, one position (first layer) */ 

if(svun<"20.0) (laylout[i][laylposn]-0.0;) 

else if(sum>20.0) {laylout[i][la/lposn]-1.0;) 

else (layloutfij[laylposn]“(double)(1.0/(1.0+exp(-sum)));) 

) /* Finished calculating four blocks, one position (first layer) */ 

laylposn-laylposn+1; 

) 

) /* end of receptive field window shifting loops (all positions) */ 

/* Propagate thru second layer */ 
fcr(bshlft-0; bshlft<4; bshift++){ 
la’'2posn-0; 

for(rowsthl-0; rowsthl<-279; rowsthl-rowsthl+31){ 

for(col8thl-rowsthl; colsthl<-rowsthl+26; colsthl++)( 
lndex-0; 

for(i-colsthl; l<-colsthl+124; l-i+31){ 
for(j-l; J<-(l+4); j++){ 
pharray[index]-laylout[bshlft][J); 
index-lndex+1; 

) 

) 

suffl-0.0; 
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for(1-0; 1<25; 1++){ 
sum-sum-i-pharray [ 1 ]; 

} 

lay2out[bshift][lay2posn] - auin*lay2vgt; 
lay2posn - lay2posn -t- 1; 

) 

1 

) 

/* Load Input array for multilayer perceptron */ 
loadcnt-0; 
for(1-0; 1<4; !++)( 
for(j-0; j<270; j++)( 

Input[loadcnt] - lay2out[1][J]; 
loadcnt-loaIcnt+l; 

) 

) 

/* Test printout module for Input array */ 

/*prlntf("flag-%d\n",flag[randnum]); 

for(1-0; KIOSO; !++){ 

printf("lnput-%f\n",Input[1]); 

) 

prlntf ("class-%d\n" .class(randntiml) ;*/ 

/* Propagate thru multilayer perceptron */ 

for(n-0; n<mldnodes; n++){ 
sum-O.0; 

for(1-0; KIOSO; i++){ 

8\im - s\un + mldwgts(n] [l]*input[l] ; 

) 

if ((sum+midthet8(n])<-20.0) {ymld[n]-0.0:) 
else If ((sxim+ffildtheta [n] )>20.0) {ymid[nl“l. 0; ) 
else {ymld[n]-(double)(1.0/(1.0+exp(-(svun+midtheta[n]>)));) 
) /* End (n) loop changing hidden layer nodes */ 

for(n-0; n<4: n-M-) { 
sum-0.0; 

for(1-0; l<mldnodes; !++){ 

sum - sum + outwgts[n][l]*ymld[1]; 

) 

If ((sum-Kmttheta[n])<-20.0) (yout(n]-0.0; ) 

else If ((siim+outtheta[n])>20.0) (yout(n]-l.0; ) 

else {yout[nl“(double)(1.0/(1.0+exp(-(sum+outtbcta[n])))),} 

) /* End (n) loop changing output layer nodes */ 

/* End propagation thru network */ 

/* Backpiopagatlon to update perceptron weights */ 

for(n“0; n<4; n++)( 

if (class[ ran dnum] — n) (dfnl-1.0;) 
else (dfnj-O.O;) 

) 
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for(n"0; n<4; n++){ 

errorout[n] « yout[n]*(l.0-yout[n])*(d[n]-yout:[n]); 

) 

for(n-0; n'Cnldnodes; n++){ 
err8\m-0.0; 

for(i-0; 1<4; !++){ 

errsum - errstun + errorout[l]*outwgtsIi][n]; 

) 

errozinld[n] ■ ymld[n]*(l.0-ymld[n])*errstun; 

) 

for(1-0; i<4; !++)( 

for(j-0; j<nldnodes; j++){ 
outwgts[i][j] - outwgts[1][j] + 

(eta/(double)mldnodes)*errorout[i]*ymid[j]; 

) 

) 

for(i-0: 1<4; !++){ 

outtheta[i] - outtheta[i] + (eta/<double)midnodes)*errorout[1]; 

) 

for (1-0; Kmldnodes; !++){ 
for(j-0; j<1080; J++){ 

mldwgts[1][J] - fflldwgts[1]{jj + (eta/1080.00)*errormid[i]*input|j); 

) 

) 

for(1-0; l<mldnodes; !++){ 

mldtheta[l] -mldtheta[l] + (eta/1080.00)*errormid[1]; 

) 

/* end of backprop */ 

) /* end of training loop (tmloop) */ 

/* Begin test loop to run thru each test vector */ 
numrlght-0; 

for(testvec-numexvec; testvec<(nuinexvec+numttvec) ; testvec-s+X 

/* Begin propagation thru network */ 

/* Begin receptive field window shifting loop */ 
laylposn-0; 

for (rowstplx-O; rowstplx<-6656; rowstpix-rowstplx <-512) { 

for(colstplx-rowstplx; colstplx<-rowstplx+120; colstpix-colstpix44)( 
lndex-0; 

for(l-col8tplx; l<-colstpix+896; i-i+128){ 
for(J-l; J<-(l+7); j++)( 

rfarray[Index]-(double)x[(testvec*8064)+J]; 

Index-lndex-fl, 

) 

) 
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/* LAMBERTIZATION AND CONTRAST NORMALIZATION MODULE (processes each 
receptive field window) */ 
absum-O.0; 

for(i-0; i<64; !++){ 
absuin-absuin+rfarray[i]; 

} 

avbright-abs\un/64.0; 
for(i-0; i<64; !++){ 
rfarray [ i]-«rfarray[i] -avbright; 

) 

esum-0.0; 

for(i-0; i<64; i++){ 
esum"»esiim+(r£array[ i]*rfarray [i] ) ; 

} 

energy»sqrt(esum); 

if (energy !- 0.0) { /* Do not normalize if energy=0.0 */ 

for(1-0; i<64; 1++){ 
rfarray!l)-rfarray[i]/energy; 

) 

) /* end normalization module */ 

/* Propagate thru i;o first hidden layer outputs */ 

for(i-0; i<4; i++)[ /* Load each 1st layer block from single window */ 

sum-0.0; 

for(j-0; J<64; j-n-){ 

sum - sum + gabor!i][j]*rfarray[j] ; 

) /* Finished calculating one block, one position (first layer) */ 

if(sum<-20.0) {laylout[i][laylposn]-0.0;) 

else if(sum>20.0) {laylout[1][laylposn]-!.0;) 

else {laylout[1][laylposn]-(double)(1.0/(1.0+exp(-sum)));) 

) /* Finished calculating four blocks, one position (first layer) */ 

laylposn-laylposn+l; 

) 

) /* end of receptive field window shifting loops (all positions ) */ 

/* Propagate thru second layer */ 
for(bshlft-0; bshlft<4; bshlft++)( 
lay2potn-0; 

for(rowsthl-0; rowsthl<-279; rowsthl-rowsthl+31){ 

for(colsthl-rowsthl; colsthl<-rowsthl+26; colsthl++)( 
lndex-0; 

for(l-colsthl; .l<-colsthl+124 ; 1-1+31) { 
for(j-i; J<-(l+4); J++)( 
pharray!index]-laylout.[bshlft] [J]; 
index-index+1; 

) 

} 

s\iffl-0.0; 
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for(l-0; i<2b; i++){ 
sum-sum+pharray[l]; 

) 

lay2out[bshift][lay2posn] - sum*lay2wgt; 
lay2posn lay2posn + 1; 

) 

) 

) 

/* Load input array for multilayer perceptron */ 
loadcnt-0 ; 
for(i«0; i<4; 1++){ 
for(j-0; j<270; j++)( 
input[loadcnt] - lay2out[1][j]; 
loadcnt-loadcnt+1; 

) 

) 


/* Propagate thru multilayer perceptron V 

for<n-0; n<liiidnodes ; n++) ( 
sujn-0.0; 

for(l-0; i<1080; 

sura ■■ s\im + midwgts[n] [i]*input [i]; 

) 

if ((sum+midtheta[n])<-20.0) {ymid(n]-0.0;) 

else if ( (sum+midtheta [n] )>20.0) 'ymid[n]>-l. 0;) 

else {ymld[n]-(double)(1.0/(1.0+exp(-(sum+midtheta[n])>));) 

} /* End (n) loop changing hidden layer nodes */ 

for(n-0; n<4; n++){ 
sum-0.0; 

for (i-0; Kmidnodes; i++) { 

sum - sum + outwgts[n][i]*ymid[l]; 

) 

if ((sum+outtheta[n])<-20.0) {yout[n]-0.0:) 

else if ((sum+outtheta[n])>20.0) {yout[n]-l.0;) 

else (yout[n]-(double)(1.0/(1.0+exp(-(sum+outtheta[n]))));) 

) /* End (n) loop changing output layer nodes */ 

/* End propagation thru network */ 

/* Decision criteria for correct classification */ 
neurtcnt-O; 
for(n“0: n<4; n++){ 

if((yout[n]>0.8) && (class[testvec]—n)) {neurtcnt-neurtcnt+1;) 
if((yout[n]<0.2) && (class[testvec]l-n)) (neurtcnt-neurtcnt+l;) 
) /* End loop for checking all four output neurons */ 
if(neurtcnt—4) (numright-numright+l;) 

) /* end of testing loop (testvec) */ 
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itnum - bigloop*numtm; 

accuracy - (float)numrlght/(float)numttvec: 

printf("Training Epoch #: %d Accuracy-%f\n", itnum, accuracy) 

) /* end of bigloop */ 

) /* end of main */ 
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f* This preprocessing progrem propagates all Image Inputs from 
the Input data set through two hard-wired layers of the 
biologically-motivated Image classification network. The 
program transforms input pixel data vectors to extracted 
feature vectors. The feature vectors, along with flag 
and class designating Integers, are written to an output 
file which can be accessed later by a separate 
classification routine. Name: FRSFROC.C */ 

#include <stdio.h> 

#include <Bath.h> 

#deflne numvec 88 /* Specify total ntunber of vectors */ 

#define numfeat 8064 /* Specify ntimber of features per vector */ 

#deflne pi 3.14159265358979 

char *malloc(); 

main(){ 

int *flag, *x, *class, gabcount; 
int 1, j, bigloop, bshlft; 

int rowstpix, colstpix, rowsthl, colsthl, index, laylposn, lay2posn 
double rfarray[64], absum, avbright, esum, energy, lay2wgt; 
double sum, fx, fy, a, b; 

double laylout[4][434], lay2o>it[4][270], pharray[25]; 
double gabor[4][64]: 

/* Open and read in data file to arrays flag[], class]], x[][] */ 
FILE •*ffp; 

fp-fopen("bigimage","r"); 
if<fp--NULL) exit(O); 

flag - (int *) malloc(numvec * sizeof(int)); 
class - (int *) malloc(numvec * sizeof(int)); 

X " (int *) malloc(numvec * numfeat * sizeof(int)); 

for(i"«0; l<n\imvec, l++)( 
fscanf (fp,"%d" ,6.f lag] 1]) ; 
for(j“0; J-Ciumfeat; J++>{ 
fscanf(fp,"%d",6ut[i*numfeat+J]); 

) 

fscanf (fp. "%d” ,6iclass[ i] ) ; 

) 

fclose(fp); 

/* Test printout module: */ 

/*for(l»«0: iCnumvec; !++)( 

printf("Flag-%d Cl8ss-%d\n",flag[i],class[l]); 

) */ 
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/* Establish constant "weight" values for Gabor function orientations 
and "phase synchronizing" summation layer */ 

lay2wgt“1.0/25,0; /* Constant weight for all 2nd hidden layer nodes >/ 

/* Initialize four Gabor function weight arrays for first layer 
(maximum value cf each function is 0.5 to aid sigmoid computation) */ 

fx-0.0; 
fy-1.0; 
gabcount-0; 
for(1-0; i<-7; !++){ 
for(J-0; J<-7; j++){ 
a- 0.2*(float)l - 0.7; 
b- 0.2*(fiost)J - 0.7; 
gabor[0][gabcount]- 

0.5*exp(-pl*(a*a/3.0+b*b/3.0))*cos(2,0*pi*(fx*a+fy‘*b)) ; 
gabcount-gabcounC+1; 

) 


fx-1.0; 
fy-0.0; 
gabcount-0; 
for(1-0; l<-7; !++){ 
for(j-0; j<-7; j++){ 
a- 0.2*(flcat)i - 0.7; 
b- 0.2*(float)J - 0.7; 
gabor[1][gabcount]- 

0.5*exp(-pi*(a*a/3.0+b*b/3.0))*cos(2.0*pl*(fx*a+fy*b)); 
gabcount-gabcount+1; 

) 


s:x-l.A142126; 
fy-i.4142136; 
gabcount-0; 
for(1-0; i<-7; 1++){ 
for(J-0; J<-7; 
a- 0.2*(float)l - 0.7; 
b- 0.2*(float)j - 0.7; 
gabor[2]Igabcount1- 

0.5*exp(-pl*(a*a/3.0+b*b/3.0))*cos(2.0*pi*(fx*a+fy*b)); 
gabcount-gabcount+1; 

) 


fx-1.4142136; 
fy-0.0-1.4142136; 
gabcount-0; 
for(l-0; i<-7; !++)( 
for(J-0; J<-7; J++){ 
a- 0.2*(float)l - 0.7; 
b- 0.2*(float)j - 0.7; 
gabor[3][gabcount]- 
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0.5’vexp( -pi*(a*a/3 • 0+b*b/3.0) )*cos (2.0*pi*(fx*a+fy*b) ) ; 
gabcount»gubcount+l; 

) 

) 

/* End of module establishing weights. Begin loop thru all vectors. */ 

for(blgloop~0; blgloop<numvec: blgloop-i-f) { 

/* Begin propagation thru network */ 

/* Begin receptive field window shifting loop */ 
laylposn-0; 

for(rowstplx-0; rowstplx<~6656; rowstpix-rovstplx+512)( 

for(col8tplx-rowstplx; colstplx<-rowstplx-«-120; colstplx-colstpix+A) { 
lndex~0; 

for<l««colstpix; l<-coistplx+896; 1-1+128){ 
for(j-l; J<-(i+7); j++)( 

rfarray[index]«(double)xl (bigloop''f8064)+j ]; 

Index-lndex+l; 

) 

) 

/* LAMBERTIZATION AND CONTRAST NORMALIZATION MODULE (processes each 
receptive field window) */ 
absvun-0.0; 

for(1-0; 1<64; !++){ 
absum-absum4rfarray[l]; 

) 

avbright-absum/64.0; 
for(i-0; i<64, !++){ 
rfarray[1]-rfarray[1]-avbrIght; 

> 

esum-0.0; 

for(i-0; i<64; !++)< 
esum-esum+(rfarray[l]*rfarrayl1)); 

} 

energy-sqrt(esuffl); 

If (energy I- 0.0) { /* Do not normalize if energy-0.0 */ 

for(i-0; 1<64: !++)( 
rfarray[1]-rfarray[1j/energy; 

) 

) /* end normalization module */ 

/* Propagate thru to first hidden layer outputs */ 

for(1-0; 1<4; l4+){ /* Load each 1st layer block from single window */ 

stun-0.0; 

for(J-0; j<64; j++){ 

sum - sum + gabor[i][J]*rfarraylJj; 

) /* Finished calculating one block, one position (first layer) */ 
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If (s\im<-20.0) (laylout[l] [laylposn]-0.0:) 

else lf(sun>20.0) (laylout[IJ[laylposn]««1.0;) 

else (laylout[i] [laylposn]-(double)(1.0/(1.0-fexp(-suin))) ;) 

) /* Finished calculating four blocks, one position (first layer) */ 

laylposn-laylposn+l; 

) 

} /* end of receptive field window shifting loops (all positions) */ 

/* Propagate thru second layer */ 
for(bshift-0; bshlft<4; bshlft++)( 
lay2posn-0; 

for(rowsthl"0; rowsthl<**279; rowsthl-rowsthl+31){ 

for(colsthl»rowsthl: colsthl<-rowsthl-t-26; colsthl++)( 

indexoO; 

for(l-colsthl; l<-colsthl+124; i-i+31)( 
for(j-l; j<-(i+4): j++){ 
pharray[lndex]-laylout[bshift][j]; 
lndex»index+l; 

) 

) 

sum-O.0; 

for(i-0; i<25; !++){ 
sum-sum+pharray[i]; 

) 

lay2out[bshift] [lay2posn] ■ suin*lay2wgt; 
lay2posn •• lay2posn + 1; 

) 

} 

) 

/* End propagation thru output of summation layer */ 

/* Print to output file: flag, new component values, 
and class (for each vector) */ 

prlntf("%d\n",flag[blgloop]); 
for(i-0; 1<4; !++){ 
for(J-0; J<270; j++){ 
prlntf("%1.16f\n",lay2out[1][J]); 

} 

) 

printf("ld\n",classIblgloop]); 


} /* end of blgloop */ 
) /* end of main */ 









f* This program inputs a data set consisting of previously 
extracted feature vectors such as those transformed by 
the biologically-motivated image processing layers in 
FREFROC.C. This routine vertically normalizes the 
input feature vectors and feeds them into a multilayer 
perceptron classifier with a variable number of hidden- 
layer nodes. Input and interim values are 
double precision. 

Name: MLFNORH.C */ 


#include <stdio.h> 
#include '<miath.h> 


#define 

numexvec 

44 

/* 

Specify 

#deflne 

numttvec 

44 

/* 

Specify 

#define 

numfeat 

1080 

/* 

Specify 

#define 

midnodes 

30 

/* 

Specify 

#define 

numruns 

50 

/* 

Specify 

#define 

nuffltm 1000 

/* 

Specify 

#define 

eta 1.0 


/* 

Specify 


of exemplar vectors */ 
of test vectors */ 
of features per vector */ 
of hidden layer nodes */ 

of result plot points */ 

training epochs between test/plot */ 

amlng rate */ 


main()( 

int *fl8g, *class, 1, j, n, neurtcnt, numrlght, itnum, randnum; 

Int blgloop, tmloop, testvec, row, col; 

double *x, outtheta[4], midtheta[midnodes], sum, errsum; 

double outwgts[4][midnodes], midwgts[midnodes][numfeat]; 

double d[4], errorout[4], errormid[midnodes], yout[4], ymid[mldnodes ] ; 

double mean, sumsqdif, s; 

float accuracy; 


/* Open and load in transformed data file of extracted features */ 
FILE *fp; 

fp-fopen("preproc.out","r"); 
if (fp—NULL) exit(O); 

flag - (int *) malloc((numexvec-^numttvec}*5izeof(int)); 
class - (int *) malloc((numexvec+numttvec)*sizeof(int)); 

X i- (double *) malloc((numexvec+numttvec)*numfeat*sizeof(double)) ; 

for(l“0; i<(nufflexvec-i-numttvec); !++){ 
fscanf(fp,"%d",&flag[l]); 
for(j«0; J<numfeat; j++)( 
fscanf (fp, "%lf" ,&x[i*numfeat-»-J ]) ; 

) 

fscanf(fp,"%d",&class[i]); 

) 

fclose(fp); 
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/* Test printout module */ 

/* for(t-0; l<(nuinexvec+numttvec); 
prlntf("%d\n",flag[l]); 
for(j-0; j<numfeat; j++){ 
printf ("il. 16f\n" ,x[l*n\jmfeat+j ]) ; 

) 

prlntf("%d\n",class[1]); 

) V 

/* Vertically normalize x[ntunexvec-t-nuinttvec] [numfeat] */ 
for (col-0; coKnumfeat; col++)( 
sum-0.0; 

for(row-0; row<nufflexvec; row++)( 
sum-sum-i'x[row*ntunfeat-fcol]; 

) 

mean-sum/(double)(numexvec); 
sumsqdlf-0.0; 

for(row-0; rou<numexvec; row++){ 

sumsqdlf-sumsqdlf+(x[ro%r*numfeat+col}-mean)*(x[row*numfeat+col)-mean) ; 

) 

s-(double)(sqrt(sumsqdlf/(double)(numexvec))); 
for(row-0; row<(nvimexvec+nuiuttvec) ; row++) { 
x[row*numfeat+col]-(x(row*numfeat+coll-mean)/s; 

) 

) /★ End column Incrementing loop */ 

/* Test printout loop */ 

/* for (row-0; row<(numexvec+n\imttvec); row++){ 
for (col-0; coKnumfeat; col++){ 

pr »tf("For row-%d col-%d normxvaiue-%f\n",row,col, 
x[row*numfeat+col]); 

> 

) */ 


srandom( (\inslgned) time (NULL)) ; /* Init randomO seed off clock */ 
/* Initialize output theta and weight arrays */ 
for(i-0; 1<4; !++)( 

outtheta[i] - (double)(random() % 101)/100.00 - 0.5; 
outtheta[i] - outtheta[i]/(double)midnodes; 

) 

for (1-0; i'Onidnodes; i++) ( 

midtheta[i] - (double)(randomO % 101)/100.00 - 0.5; 
midtheta[i] - midtheta[l]/(double}numfeat; 

) 

for(i-0; 1<4; 1++)( 

for(J-0; J<3aidnodes; J++){ 

outwgts[i][j] - (double)(randomO % 101)/100.00 - 0.5; 
outwgtsii][J] - outwgts[i)[J]/(double)midnodes; 

) 

) 
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for(i-0; l^Ooldnodes; 1++) ( 
for(J"0; J<numfeat; j++){ 

Didvgts[l][j] - (double)(randooO % I01)/100.00 - 0.5; 
mldwgtsil][J] " mldwgts[l][J]/(double)numfeat: 

I 

} 

/* End of initialisation; begin loops for training and testing. */ 
for(bigloop«l; bigloop<-numzuns; bigloop4-+){ 

/* Begin training loop */ 
for(tmloop»0; tmloop<n\jmtm; tmloop++){ 
randnum-(randoic() % numexvec); 

/* Propagate thru Oiultllayer perceptron */ 
for(n“0; nonidnodes; n++){ 
su]nx>0.0; 

for(i-0; l<numfeat; !++){ 

stun ■■ sum + midwgts[n] (l]*x[randnum*numfeat+i} ; 

) 

if ((sum+mldtheta[n])<-20.0) {ymld[n]“0.0;) 

else if ((s\un+nldtheta[n])>20.0) {ymid[n]-l.0;) 

else (ymid[n]-(doublfc)(1.0/(1.0+exp(-(sum+midtheta[n]))));) 

) /* End (n) loop changing hidden layer nodes */ 

for(n-0; n<4; n-H-){ 
suu«0.0; 

for(i"0; Kmldnodes; !++){ 

sum “ sum + outwgtsfn] ti]*yTnld[i] ; 

) 

if ((sum+outtheta[n])<-20.0) {yout[n]-0.0;) 

else if ((s’jua+outthetalnl)>20.0) {yout[nj«1.0;) 

else {yout[n]‘«(double)(1.0/(1.0+exp(-(sum+outtheta[n ]))));) 

} /* End (n) loop changing output layer nodes */ 

/* End propagation thru network */ 

/* Bsckpropagatlon to update perceptron weights */ 

for(n~0; n<4; n++){ 

if (class[randnum] <-» n) (d[n]-1.0;} 

else (d[n]»0.0;) 

) 

for(n-0; n<4; n++)( 

errorout[n] -• yout[n]*(1.0-yout[n])*(d[n]-yout[n]); 

) 

for(n-0; n<anldnodes; n-H-){ 
errsum-O.0; 

for(i-0; i<4; !++){ 

erraum errsum + errorout[i]*outwgts[1][n]; 

) 

errormld[n] » ymid[n]*(1.0-ymidln])*errsum; 

J 
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for<1-0; i<4; !++)( 

for(j»>0; j<midnodes; j++){ 
outwgtsti] [j] - outwgts[l] [ j] 

<et«t/(double)inldnodes)*errorout[l]*ymld[ j ] ; 

) 

) 

for<1-0; i<4; !++)( 

outtheta[l] - outtheea[l] + (eta/(double)i&ldnodes)*errorout[l]; 

) 

for(1-0; l<3nidnodes; i++){ 
for(j-0; Joiumfeat; j-J’+X 
inidwgts[i] [J] - mldwgts[l] [ j] + 

(eta/(double)nu«fftat)*errorTnldtl]*x[fandnuin*ruinfeat+j ] 

) 

) 

for(1-0; Kmldnodes; l4+){ 

inldtheta[l] » mldtheta(l] + (eta/<double)nuinfeat)*errormld[ i] ; 

} 

/* end of backprop */ 

) /* end of training loop (tmloop) */ 

/★ Begin test loop to run thru each test vector */ 
numrlght-O; 

for(testvec«nuinexv9c; testvec<(nujne3cvec+n'uinttvec) ; testvec++) { 

/* Propagate thru multilayer perceptron / 

for(n“0; n<Jnldnodes; n++){ 
sum-0.0; 

for(1-0; Knumfeat; i++){ 

sum " sum + midwgtsfn][i]*x[testvec*numfeat+l]; 

) 

if ((sum+mldthetaln] )<-20.O') {yiiild[«]-0.0;) 

else if ((sum+mldcheta(n})>20.0) {yaldtn]-l.0;) 

else {ymld[n)-(dov’ble)(1.0/(1.0+exp(-(s\im+midtheta[n])))) ; ) 

) /* End (n) loop changing hidden layer nodes */ 

for(n-0; n<4; n++)( 

5um-0.0; 

for(l-0; l^ldnodes; !++){ 

Bum - sum + outwgtstn] [l]*)TBld{l]; 

J 

If ((sum+outtheta[n])<-20.0) {yout[nj-0.0;) 

else If ((8um+outtheta[n])>20.0) (youtln]-1.0;) 

else lyout[n]-(double) (1.0/(1.0+exp(-(sum+outtheti? [n] ))));) 

) /* End (n) loop changing output layer nodes */ 

/* End propagation thru network */ 
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/* Decision criteria for correct classification */ 
neurtcnt»0; 
for(n-0; n<4; n++){ 

If ((yout[n]|>0.8) && (class(testvec]—n)) (neurtcnt-neurtcnt+l;) 
lf((yout[n]<0.2) && (class[testvec]f-m)) (neurtcnt-neurtcnt+l;) 

) /* End loop for checking all four output neurons */ 

If(neurtcnt—4) (numrlght-numrlght+l;) 

) /* end of testing loop (testvec) */ 

Itniun - blgloop*nuintm; 

accuracy - (float)ntimrlght/(float)n\jmttvec; 

prlntf("Training Epoch #: %d Accuracy-%f\n", Itnum, accuracy); 
) /* end of blgloop */ 

} /* end of main */ 
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/* Image classification network version for processing input 
from PREPROC.C. Classifier incorporates a single-layer, 
second-order neural network algorithm with automatic 
vertical normalization of inputs and multiplied combos. 

Multiplied combos are of same location node outputs 
among the four different orientation groups. */ 

/* Name; IMGN0RM2.C */ 

#include <stdio.h> 

#include <math.h> 

#define numneur 4 

#define numexvec 44 
#define numttvec 44 
^define numfeat 1080 

#define numruns 500 
#define numtrn 100 

#define eta 1.0 

char *malloc(); 

main(){ 

/* flag[] holds ID number of each data vector, class[] holds training 
class associated with each data vector, x[][] holds feature 
components of data vectors, quadin[][] holds second-order multiplied 
combos of X values, wlin[][] and wquad[3[] hold weights for all inputs 
(including higher order combos), theta[] holds threshold values for 
each neuron */ 

int *flag, *class; 

int i, j, n, randnum, numquad, qcount; 
int row, col, spot; 

int neurtcnt, numright, bigloop, tcnloop, testvec, itnum; 
double d, y, error, neweta; 
double *x, *wquad, *quadin; 

double theta[numaeur], wlin[numneur*numfeat]; 
double sum, linsum, quadsum, mean, s, sumsqdif, fanin; 
float accuracy; 

/* This module opens an input file and reads data into the 
arrays flag[], class[], and x[][3 */ 

FILE *fp; 

fp=fopen{"preproc.out",”r"); /* Insert correct data file name here */ 

if(fp==NULL) exit(0); 

/* Dynamic alloc of arrays prior to loading in data file values */ 

flag = (int *) malloc((numexvec+nurottvec) * sizeof(int)); 
class = (int *) malloc((numexvec+numttvec) * sizeof(int)); 

X = (double *) malloc((numexvec+numttvec)*numfeat * sizeof(double)); 


/* Specify total number of neurons across slab */ 
/* Number of neurons = number of classes */ 

/* Specify number of exemplar vectors */ 

/* Specify number of test vectors */ 

/* Specify number of features per vector */ 

/* numfeat also represents # FEATURE COLUMNS */ 

/* Specify number of result plot points */ 

/* Specify # training iters between test/plot */ 

/* Specify training factor */ 

/* Must define for use of malloc thru program */ 
/*Need for this statement is compiler-dependent*/ 
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for(i=0; i<(numexvec+numttvec); ++i){ 
fscanf (fp,"%d'',&flag[i]) ; 

for(j=0; j<nuinfeat; ++j){ 
fscanf (fp, '•%lf ", &x[ i*nuinfeat+j ]); 

} /* end inner loop */ 

fscanf (fp, ''%d'', &class[i]) ; 

} /* end outer loop */ 

fclose(fp); 

/* This test module prints out input data containing first two and 
final feature components of each vector; checks the above file 
opening module (comment out for final program) */ 

/* for(i=0; i<(numexvec+nuffittvec); ++i){ 

printf("For flag=%d, class=%f, check features are %f %f %f\n", 
flag[i], class[i], x[i*numfeat+0], x[i*numfeat+l], 
x[i*numfeat+(numfeat-1)]); 

) */ 

/* This statement establishes the number of second-order 

combos of feature inputs FOR EACH VECTOR (without commutative 
redundancy), unique to this problem. Will be used for array 
allocation. */ 

numquad=2700; 

/* Fill 2-D arrays for second order input combos; 
quadin[numexvec+numttvec][numquad] */ 

quadin=(double *) malloc((numexvec+numttvec)*numquad*sizeof(double)) 
for(row=0; row<(numexvec+numttvec); row++){ 
qcount=0; 

for(spot=0; spot<270; spot++){ 

for(i=spot; i<=(spot+810); i=i+270){ 
for(j=i; j<=(spot+810); j=j+270){ 

quadin[row*numquad+qcount] = x[row*numfeat+i] * x[row*numfeat+j] 
qcount = qcount +1; 

) 

) 

) 

} /* end vector row incrementing loop */ 

/* Test printout module */ 

/* for(row=0; row< (numexvec+numttvec) ; row-t r) { 
for(col=0; col<numquad; col++){ 

printfC'For row=%d col=%d quadin=%f\n",row,col, 
quadin[row*numquad+col]); 

) 

) */ 
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/* Begin VERTICAL NORMALIZATION of x, quadin arrays */ 


/* Vertically normalize x[numexvec+numttvec][numfeat] */ 
for(col=0; col<numfeat; col++){ 
sum=0.0; 

for(row=0; row<numexvec; row++){ 
sum=sum+x[row*numfeat+col]; 

) 

mean=sum/(double)(numexvec); 
sumsqdif=0.0; 

for(row=0; row<numexvec; row++){ 

sumsqdif=sumsqdif+(x[row*numfeat+col]-mean)*(x[row*numfeat+col]-mean) 

} 

s=(double)(sqrt(sumsqdif/(double)(numexvec))); 
for(row=0; row<(numexvec+numttvec); row++){ 

X[row*numfeat+col]=(x[row*numfeat+col]-mean)/s; 

} 

) /* End column incrementing loop */ 

/* Test printout loop */ 

/* for(row=0; row<(numexvec+numttvec); row++){ 
for(col=0; col<numfeat; col++){ 

printfC'For row=%d col=%d normxvalue=%f\n",row,col, 

X[row*numfeat+col]); 

) 

) */ 

/* Vertically normalize quadin[numexvec+numttvec][numquad] *, 
for(col=0; col<numquad; col++){ 
sum=0.0; 

for(row*0; row<numexvec; row++){ 
sum=sum+quadin[row*numquad+col]; 

) 

mean=sum/(double)(numexvec); 


sumsqdif=0.0; 

for(row=0; row<numexvec; row++){ 
sumsqdif=sumsqdif+ 

(quadin[row*numquad+col]-mean)*(quadin[row*numquad+col]-mean); 

) 

s=(double)(sqrt(sumsqdif/(double)(numexvec))); 
for(row=0; row<(numexvec+numttvec); row++){ 
quadin[row*numquad+col]=(quadin[row*numquad+col]-mean)/s; 

) 

) /* End column incrementing loop */ 

/* Test printout Ir-op */ 

/* for(row=0; row<(numexvec+numttvec); row++){ 
for(col=0; col<numquad; col++){ 

printfC'For row=%d col=%d normquadvalue=%f\n",row,cc , 
quadin[row*numquad+col]); 

, i/ 

/* End of vertical normalization */ 
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srandoin( (unsigned) time (NULL)) ; /* Init random() seed off clock */ 

/* Loop to initialize theta array (one for each neuron) with random 
floating point value between -0.5 and 0.5 */ 

for(i=0; i<numnGar; i++){ 

theta[i]=(double)(random() % 101)/1C0.00 - 0.5; 

/* printfC'For i=%d theta=%f\n'', i, theta[i]) ; */ 

} 

/* Loop to initialize linear input weights array wlin[numneur][numfeat] 
with random floating point values between -0.5 and 0,5 */ 

for(n=0; n<numneurr n++)( 

/* printfC'For neuron # %d\n",n); */ 
for(i-0; i<numfeat; i++){ 

wlin[n*numfeat+i]=(double)(randomO % 101)/100.00 - 0.5; 

/* printfC'For i=%d wlin=%f\n'', i,wlin[n*numfeat+i]) ; */ 

} 

} /* end of nested loop */ 

fanin = (double)numfeat + (double)numquad; 

/* fanin = 1; */ 

neweta == (double) eta / fanin; 

/* printf("neweta= %f\n",neweta); */ 

/* Loop to initialize quadratic input weights array 

wquad[numneur][numquad] with random floating point values between 
-0.5 and 0.5 */ 

wquad = (double *) malloc(numneur*numquad*sizeof(double)); 
for(n=0; n<numneur; n++){ 

/* printfC'For neuron # %d\n",n); */ 
for(i--0; i<numquad; x++) { 

wquad[n*numquac? + i] = (double) (randOin() % 101)/100.00 - 0.5; 

/* printfC'For i=id quadweight=%f\n",i,wquad[n*numquad+i]); */ 

) 

) /* end of nested loop */ 

/* End of initialization; Begin loops for training and testing */ 
for(bigloop=l; bigloop<=numruns; bigloop++){ 

/* Begin training loop */ 

for (trnloop“0; trnloop<numtrn; trnloop-f + ) { 

/* Randomly select an exemplar vector row (flag #) */ 
randnum=(random() % numexvec); 

/* printf("Randomly selected vector flag number is %d\n",randnum); */ 

for(n-=0; n<numneur; n++) { /* Loop to update weights for each neuron */ 

linsum=0.0; 

for(i*=0; i<numfeat; i++) ' 

linsum=linsum+wl in [ n*numf eat+i ] *x [ randnum*numf eat-i i ] ; 

} 

quadsum=0. 0 ; 
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for(i=0; i<numquad; i++){ 

quadsum=quadsuin+wquad [ n*rur!\quad+i ] *quadin [ randnuin*nuitiquad+ i ] ; 

} 

suin=linsuin+quadsuin; 

if ((SU 1 I 1 +theta [n]) <-20.0) {y=0.0;} 

else if ((suin+theta[n]) >20.0) {y=1.0;) 

else {y= (double) (1.0/ (1.0+exp(-(suin+theta[n]))));) 

/* Criteria for training each class to fire a specific neuron */ 
if (class[randnuin]==n) {d=1.0;) 
else {d=0.0;} 

error=y*(1.0-y)*(d-y); 

/* printf("error=%f\n",error); */ 
for(i=0; i<numfeat; i++){ 

wl in [ n*nuinf eat+i ] =wl in [ n*nuinf eat+i ] +neweta*error *x [ randnum*nuiT\f eat+ i ] 
/* printfC'For i=%d wlin=%f\n'', i,wlin[n*nuinf eat+i]) ; */ 

} 

for(i=0; i<nuinquad; i++) { 

wquad [ n*numquad+i ] =wquad [ n*nuinquad+i ] + 

neweta*error*quadin[randnuin*nuinquad+i ] ; 

/* printfC'For i=%d quadweight=%f\n", i, wquad [ n*nurnquad+i ]) ; */ 

} 

theta[n]=theta[n]+neweta*error; 

/* printfC'For neuron # %d theta+%f\n”,n,theta[n]); */ 

) /* end of weight training loop for each neuron's weights */ 

) /* end of complete training loop */ 

/* Begin test loop to run thru each test vector */ 
nuinright=0; 

for(testvec=^nuinexvec; testvec<(numexvec+numttvec) ; testvec++){ 

/* printf("testvec=%d\n",testvec); */ 

neurtcnt=0; 

tor(n=0; n<numneur; n++){ 
linsuin=0.0; 

for(i=0; i<nuinf6it; i++) { 

linsuin=linsum+wlin [ n*nuinfeat+i ] +x [ testvec*nuinfeat+ i ] ; 

) 

quadsuin=0.0; 

for(i=0; i<nuinquad; i++) { 

quadsuIn=quadsum-^wquad [n*nui!iquad+i ] *quadin [ testvec*numquad+ i ] ; 

} 

sum=l insuiti+quadsum; 

if ((sum+theta[n])<-20.0) (y=0.0;) 

else if ((sum+theta[n])>20.0) (y=1.0;> 

else {y=(double) (1.0/ (1.0+exp(-(suiii+theta[n] )))),*) 

/* printfC*y=%f class[testvec]=%f\n",y,class[testvec]); */ 
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/* Test decision criteria for identifying classes */ 
if ((y>0.8) && (class[testvec]==n)) ineiirtcnt=neurtcnt+l?) 
if ((y<0.2) (class[testvec] i^n)) {neurtcnt=neurtcnt+l;} 

) /* end of neurons (n) correct test loop */ 

if (neurtcnt==nuinneur) {numrighV.=miSiright+l;} 

) /* end of testing loop running thru each test vector (testvec) */ 

itnum=bigloop*nuintrn; 

accuracy- (float) numright/ (f loat) nuinttvec; 

printf("Training iteration #: %d Accuracy=%f\n",itnum,accuracy); 

) /* end of bigloop */ 

) /* end of main */ 
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